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PREFACE 


This  report  is  published  to  provide  coastal  engineers  with  a  mathematical 
modeling  procedure  for  predicting  shoreline  evolution  resulting  from  the  con¬ 
struction  of  navigation  and  shore  structures.  The  model  is  calibrated  with  a 
test  case  at  Holland  Harbor,  Michigan.  This  report  is  a  continuation  of  an 
investigation  by  Le  Mehaute  and  Soldate  (1977)  to  determine  the  feasibility 
of  applying  numerical  models  to  real  case  situations.  The  work  was  carried 
out  under  the  coastal  structures  program  of  the  Coastal  Engineering  Research 
Center  (CERC) . 
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CONVERSION  FACTORS,  U.S.  CUSTOMARY  TO  METRIC  (SI)  UNITS  OF  MEASUREMENT 


U.S.  customary  units  of  measurement  used  in  this  report  can  be  converted  to 
metric  (SI)  units  as  follows: 


Multiply 

by 

To  obtain 

inches 

25.4 

millimeters 

2.54 

centimeters 

square  inches 

6.452 

square  centimeters 

cubic  inches 

16.39 

cubic  centimeters 

feet 

30.48 

centimeters 

0.3048 

meters 

square  feet 

0.0929 

square  meters 

cubic  feet 

0.0283 

cubic  meters 

yards 

0.9144 

meters 

square  yards 

0.836 

square  meters 

cubic  yards 

0.7646 

cubic  meters 

miles 

1.6093 

kilometers 

square  miles 

259.0 

hectares 

knots 

1.852 

kilometers  per  hour 

acres 

0.4047 

hectares 

foot-pounds 

1.3558 

newton  meters 

millibars 

1.0197  x  10“3 

kilograms  per  square  centimeter 

ounces 

28.35 

grams 

pounds 

453.6 

grams 

0.4536 

kilograms 

ton,  long 

1.0160 

metric  tons 

ton,  short 

0.9072 

metric  tons 

degrees  (angle) 

0.01745 

radians 

Fahrenheit  degrees 

5/9 

Celsius  degrees  or  Kelvins1 

JTo  obtain  Celsius 

(C) 

temperature  readings 

from  Fahrenheit  (F)  readings. 

use  formula:  C  = 

(5/9)  (F  -32). 

To  obtain  Kelvin 

(K) 

readings,  use  formula 

K  =  (5/9)  (F  -32)  +  273.15. 
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SYMBOLS  AND  DEFINITIONS 


a 

B 

C£ 

Co 

d 

dVl 

dV2 

% 

% 

g 

Hb 
H o 

kd 

kr 

h 

k 

Lb 
^ o 

l 

M 

m 

ox 

oy 

02 

Q 

% 


proportionality  constant  (6  =  aa) 

height  of  bluff  (in  case  of  erosion)  or  berm  (in  case  of  accretion) 

phase  velocity  at  breaking 

deepwater  phase  velocity 

water  depth  at  yQ 

water  depth 

breaking  depth 

variation  of  sand  volume  due  to  horizontal  displacement  during  time,  dt 

variation  of  sand  volume  due  to  beach  slope  variation  during  time,  dt 

group  velocity  at  breaking 

deepwater  wave  group  velocity 

gravity  acceleration 

breaking  wave  height 

deepwater  wave  height 

diffraction  coefficient 

refraction  coefficient  from  deep  water  to  the  line  of  breaking  inception 

percentage  of  silt 

time-constant  parameter 

wavelength  at  breaking 

deepwater  wavelength 

length  of  groin  perpendicular  to  shore 
loss  of  sand  by  rip  currents 
beach  slope  at  shoreline 

horizontal  axis  parallel  to  the  average  initial  beach  profile 
horizontal  axis  perpendicular  to  x,  oriented  positively  seaward 
vertical  axis,  positive  upward 
longshore  energy  flux 

dimensionless  littoral  drift  discharge 
littoral  drift  discharge  in  cubic  yards  per  year 
loss  of  bluff  volume  (in  case  of  erosion)  due  to  silt 
loss  of  bluff  volume  (in  case  of  erosion)  dee  to  silt 
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SYMBOLS  AND  DEFINITIONS--Continued 
loss  of  sand  by  wind 

sand  discharge  in  a  direction  perpendicular  to  the  beach 
average  bottom  slope  to  depth  of  no  motion 

horizontal  displacement  at  depth,  Dc,  characterizing  slope  variation 

shoreline 

wave  period 

time 

volume  of  sand  over  a  stretch  of  shoreline.  Ax,  unity 
shoreline  coordinates 
dimensionless  shoreline 
deepwater  limit  of  the  beach 

function  of  x,  y,  t  defining  the  bottom  topography 
shoreline  elevation  above  a  datum  z  =  o 

7 . 5  x  lo3  fl  x  T 
64 

deepwater  wave  angle  with  the  ox-axis 
angle  of  breaking  with  the  shoreline 
deepwater  wave  angle  with  the  shoreline 

angle  of  the  beach  limit  with  the  ox-axis  in  front  of  a  groin 

horizontal  displacement  of  the  beach  profile  during  time,  dt 

function  characterizing  beach-slope  variation  as  function  of  x 

function  characterizing  beach-slope  variation  as  function  of  t 

function  characterizing  beach-slope  variation  with  time  due  to  rapid 
change  of  sea  level 

function  characterizing  beach-slope  variation  with  time  due  to  manmade 
structures  such  as  groin  or  navigation  works 

a  symptotic  (in  axium)  value  of  ,\y"  when  t  >  ■» 

variation  of  wave  direction  by  diffraction 

angle  of  wave  ray  with  a  perpendicular  to  shore  in  a  diffraction  zone 
water  density 
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A  NUMERICAL  MODEL  FOR  PREDICTING  SHORELINE  CHANGES 


by 

Bernard.  Le  Mehaute  and  Mills  Soldate 

I.  INTRODUCTION 

This  report  establishes  a  mathematical  model  for  shoreline  evolution  and 
calibrates  the  model  with  a  test  case  at  Holland  Harbor,  Michigan.  Even  though 
the  mathematical  model  is  general  and  could  be  applied  to  a  number  of  situa¬ 
tions,  the  emphasis  is  on  the  Great  Lakes,  and  more  specifically,  on  shoreline 
evolution  near  navigation  structures. 

An  interim  report  by  Le  Mehaute  and  Soldate  (1977)  reported  on  the  feasi¬ 
bility  of  applying  existing  mathematical  models  to  real  case  situations.  This 
report  is  a  continuation  of  that  first  investigation. 

The  present  mathematical  model  includes  many  of  the  characteristics  already 
covered  in  the  literature.  In  addition,  the  model  presents  an  integrated  ap¬ 
proach  on  a  large  number  of  phenomena  previously  neglected.  Its  main  purpose 
is  to  develop  a  practical  numerical  scheme  which  could  be  used  to  predict  shore¬ 
line  evolution.  The  model  would  then  be  able  to  point  out  the  shortcomings  in 
the  present  state  of  knowledge.  Therefore,  the  mathematical  model  covers  some 
aspects  of  shoreline  evolution  which  cannot  be  quantified  with  the  data  obtained 
from  the  considered  test  case  of  Holland  Harbor.  The  mathematical  model  then 
has  to  be  regarded  as  a  research  guide  for  the  future. 

One  important  aspect  is  the  effect  of  sand  size  and  density.  It  is  well 
known  that  the  rate  of  shoreline  erosion  and  sediment  loss  is  largely  affected 
by  these  parameters.  The  fine  sand  is  transported  offshore,  while  larger  size 
sand  tends  to  proceed  alongshore  according  to  a  littoral  drift  formula.  This 
effect  could  not  presently  be  quantified;  therefore,  it  is  introduced  in  the 
mathematical  model  as  a  constant. 

The  present  mathematical  model  can  continuously  be  upgraded  as  the  state- 
of-the-art  progresses,  and  as  the  model  is  tested  for  a  large  number  of  cases. 

It  is  important  to  remember  that  the  model  deals  with  long-term  shoreline  evo¬ 
lution  as  defined  in  Le  Mehaute  and  Soldate  (1977).  Short-term  evolution  must 
be  considered  as  local  perturbations  which  are  superimposed  on  the  presently 
defined  topography. 

Three  time  scales  of  shoreline  evolution  which  can  be  distinguished  are 
(a)  geological  evolution  over  hundreds  and  thousands  of  years,  (b)  long-term 
evolution  from  year-to-year  or  decade,  and  (c)  short-term  or  seasonal  evolu¬ 
tion  during  a  major  storm. 

Associated  with  these  time  scales  are  distances  or  ranges  of  influence  over 
which  changes  occur.  The  geological  time  scale  deals,  for  instance,  with  the 
entire  area  of  the  Great  Lakes.  The  long-term  evolution  deals  with  a  more 
limited  stretcli  of  shoreline  and  range  of  influence;  e.g.,  between  two  head¬ 
lands  or  between  two  harbor  entrances.  The  short-term  evolution  deals  with 
the  intricacies  of  the  surf  zone  circulation;  e.g.,  summer-winter  profiles, 
bar,  rhythmic  beach  patterns,  etc. 
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For  the  problem  under  consideration,  long-term  evolution  is  of  primary 
importance,  the  short-term  evolution  appearing  as  a  superimposed  perturbation 
on  the  general  beach  profile.  Evolution  of  the  coastline  is  characterized  by 
low  monotone  variations  or  trends  on  which  are  superimposed  short  bursts  of 
rapid  development  associated  with  storms. 

The  primary  cause  of  long-term  evolution  is  water  waves  or  wave-generated 
currents.  Three  phenomena  intervene  in  the  action  which  waves  have  on  shore¬ 
line  evolution: 

(a)  Erosion  of  beach  material  by  short -period  seas  versus  accre¬ 
tion  by  longer  period  swells; 

(b)  effect  of  water  level  changes  on  erosion;  and 

(c)  effect  of  breakwaters,  groins,  and  other  structures. 

Although  mathematical  modeling  of  shoreline  evolution  has  inspired  some 
research,  it  has  received  only  limited  attention  from  practicing  engineers. 

The  present  methodology  is  based  mainly  on  (a)  the  local  experience  of  engi¬ 
neers  who  have  a  knowledge  of  their  sectors,  understand  littoral  processes, 
and  have  an  inherent  intuition  of  what  should  happen;  and  (b)  movable-bed 
scale  models  that  require  extensive  field  data  for  their  calibration. 

In  the  past,  theorists  have  been  dealing  with  idealized  situations,  rarely 
encountered  in  engineering  practice.  Mathematical  modelers  apparently  have 
long  been  discouraged  by  the  inherent  complexity  of  the  phenomena  encountered 
in  coastal  morphology.  The  lack  of  well-accepted  laws  of  sediment  transport, 
offshore-onshore  movement,  and  poor  wave  climate  statistics  have  made  the 
task  of  calibrating  mathematical  models  very  difficult.  Considering  the  im¬ 
portance  of  determining  the  effect  of  construction  of  long  groins  and  naviga¬ 
tion  structures,  and  the  progress  made  in  determining  wave  climate  and  littoral 
drift,  a  mathematical  approach  now  appears  feasible. 

The  complexity  of  beach  phenomena  could,  to  a  large  extent,  be  taken  into 
account  by  a  numerical  mathematical  scheme  (instead  of  closed-form  solutions), 
dividing  space  and  time  intervals  into  small  elements  in  which  the  inherent 
complexity  of  the  morphology  could  be  taken  into  account.  Furthermore,  better 
knowledge  of  the  wave  climate,  a  necessary  input,  will  allow  a  better  calibra¬ 
tion  of  coastal  constants  (such  as  those  in  the  littoral  drift  formula). 

In  past  investigations  (Le  Mehaute  and  Soldate,  1977),  many  important  ef¬ 
fects  have  been  neglected,  such  as  combined  effects  of  wave  diffraction  around 
littoral  obstacles,  change  of  sea  level,  height  of  berm  and  bluff,  beach  slope, 
etc.  The  present  investigation  attempts  to  include  all  the  important  factors 
associaied  with  long-term  shoreline  evolution.  In  the  case  of  the  Great  Lakes, 
importance  must  be  given  to  variations  in  lake  level.  The  coastal  zone  is  de¬ 
fined  by  a  three-dimensional  bottom  topography  instead  of  a  two-dimensional 
shoreline  (or  two  lines)  as  in  the  previous  cases. 

Mathematical  modeling  is  essentially  approached  by: 

(a)  Understanding  the  phenomenology  of  shoreline  evolution 

quantitatively. 
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(b)  Calibrating,  determining  empirical  constants,  and  separating 
various  effects,  such  as  due  to  change  of  lake  level,  navigation 
structures,  etc. 

(c)  Predicting  long-term  future  evolution. 

(d)  Assessing  the  effect  of  future  construction. 

The  mathematical  model  presented  in  this  report  is  just  a  first  step  in  this 
direction.  Much  of  the  model  may  be  modified  after  the  application  of  the  mathe¬ 
matical  model  to  other  cases  which  could  be  used  for  further  calibration. 

Theoretical  developments  for  the  model  are  presented  in  Section  II.  Section 
III  describes  the  shoreline  evolution  recorded  at  Holland  Harbor,  Michigan,  and 
compares  the  mathematical  model  with  the  test  case  investigated.  Section  IV 
provides  recommendations  and  conclusions.  An  Appendix  describes  a  computer 
program  to  investigate  shoreline  behavior. 

II.  THEORETICAL  DEVELOPMENTS 

Consideration  is  given  to  a  coastal  zone  limited  by  boundaries  at  a  small 
distance  from  the  surf  zone  (Fig.  1).  The  bottom  topography  is  defined  in  a 
three-coordinate  system,  oxyz,  by  a  function  Z£  =  f(x,y,t)  where  the  ox-axis 
is  parallel  to  the  average  shoreline  direction,  the  oy-axis  is  perpendicular 
seaward,  and  oz  is  positive  upward  from  a  fixed  horizontal  datum.  The  angle 
of  the  shoreline  with  the  ox-axis  is  small.  The  shoreline  is  defined  by  y  =  ys, 
z  =  zs  =  zi(x,ys,t)  which  also  defines  the  sea  level  as  a  function  of  time. 


Figure  1.  Coastal  zone. 


The  deepwater  limit  of  the  beach  is  /  =  yQ.  (This  limit  defines  the  con¬ 
tour  line  where  the  sand  is  no  longer  moved  by  wave  action.)  The  water  depth 
at  y  =  yQ  is  Da.  It  will  be  assumed  that  De  remains  constant  as  sea  level 
and  beach  profiles  change.  Therefore,  3zs/3t  =  Sz^/St. 

B  is  the  height  of  the  bluff  in  case  of  erosion,  i.e.,  when  3ys/3t  <  0, 
and  the  height  of  the  berm  in  case  of  accretion,  i.e.,  when  3ys/3t  >  0  (Fig.  2). 

EROSION  ACCRETION 


Figure  2.  Height  of  bluff  or  berm. 

The  quantity  of  sand  over  a  stretch  of  shoreline.  Ax  =  unity  and  bounded  by 
the  datum,  z  =  0,  y  =  0,  and  the  beach  profile  z^  at  time,  t,  is 


V(t) 


zfc(x,y,t)  dy. 


Assume  that,  for  some  reason,  the  beach  profile  changes  during  an  infini¬ 
tesimal  amount  of  time,  dt.  A  further  assumption  is  that  the  initial  beach 
profile  which  is  considered  at  time,  t  =  t^,  could  be  the  normal  "equilibrium 
profile."  (The  equilibrium  profile  may  never  exist  under  varying  prototype 
conditions  (similarly  two-dimensional  waves  never  exist),  but  it  is  a  conven¬ 
ient  idealized  concept  which  could  be  approached  in  two-dimensional  wave  tank 
experiments.  In  this  case,  it  could  be  defined  as  the  statistical  long-term 
average  beach  profile  which  exists  under  a  given  wave  climate.  The  model  here 
is  actually  independent  from  this  definition.) 

The  departure  and  modification  from  this  initial  beach  profile  can  be 
characterized  by  (see  Fig.  3): 

(a)  A  translation  in  the  yz  plane  4'<:ined  by  an  elementary 
vector  of  components. 


9/s  dD 
3t  ’  dt 

where  dD/dt  is  the  rate  of  change  of  sea  level.  Note  that  this 
translation  is  independent  from  the  beach  profile  and,  in  particular, 
if  the  beach  profile  normally  exhibits  a  number  of  significant  bar 
formations,  under  normal  conditions  the  translation  will  reproduce 
this  characteristic  at  the  same  water  depth. 

(I))  A  perturbation  characterizing  the  departure  or  variation 
from  the  initial  profile.  Since  the  rate  of  the  vertical  component 
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z 


Figure  3.  Change  of  beach  profile. 

of  translation  is  dD/dt,  the  perturbation  can  be  defined  only  by  a 
horizontal  displacement.  At  the  deepwater  limit,  Dc,  the  horizontal 
displacement  characteri zing  the  rate  of  change  of  the  average  beach 
slope  is  defined  by 


3S^(x,t) 

— 

so  that 

a  ^y's  ^y 

It  IF-  +  IF- 


The  quantity  of  sand  within  the  considered  domain  at  time  (t  +  dt)  is 


V(t  +  ut) 
and  the  variation 


z£>(x»y>t  +  dt)  dv 


dV 

dt 


dzy 

dy  (y1  and  0  are  fixed  limits) 


dv=  fYl  (dJb  +  +  dx\d 

dt  L  \3t  3y  dt  3x  dt  /  y 


(1) 


where  3z^/3x  is  the  variation  of  beach  elevation  along  the  ox-axis.  Since 
the  angle  of  the  beach  with  the  ox-axis  is  small 
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3z^  /3zfc 
By  ;  \3y 


m  is  the  beach  slope) 


i.e.,  the  change  in  z  occurring  along  the  beach  is  small  relative  to  change 
across  profile  and  the  beach  profile  variation  over  a  distance  dx  is  small 
but  remains  an  infinitesimal  of  higher  order  than  over  a  distance  dy. 

For  the  same  reason,  the  velocity  of  beach  variation  along  the  ox-axis 
dx/dt  is  small  as  compared  to  the  beach  variation  along  oy;  therefore 

3z£>  dx  <<;  3z&  dy 
3x  dt  <  3y~  dt 

and  could  be  neglected.  To  evaluate  the  other  terms  it  is  simpler  to  refer  to 
Figure  3.  It  is  seen  that  the  integral  is  equal  to  the  sum  of: 

(a)  The  change  of  volume  of  sand,  dv/dt,  due  to  the  translation, 
3ys/3t,  is  the  sum  of  the  change  of  volume  due  to  the  horizontal 
translation  (Fig.  4) 


(B  +  C) 


and  the  change  of  volume  due  to  the  vertical  displacement 

,  ,  dD 

dt 

where  y^  is  the  ordinate  characterizing  the  location  of  the  bluff. 


Figure  4.  Changes  of  sand  volume  due  to  translation 
of  one  beach  profile. 


14 


The  sea  level  displacement,  dD,  being  small,  the  difference  in 
the  quantity  of  sand  represented  by  the  two  triangles  ABC  and  A'B'C' 
in  Figure  4  is  considered  as  a  difference  of  infinitesimals  and  there¬ 
fore  is  neglected.  So,  combining  equations  (2)  and  (3) 


dV 

dt 


3ys 

■  (B  +  V  ar  -  (y*  - 


dD 

dt 


(4) 


It  is  seen  that  when  dV/dt  =  o 


dy,? 


yg  -  y2> 

B  +  D 


where  S  is  the  average  bottom  slope.  Note  that  this  expression  is 
independent  of  the  beach  profile  and,  therefore,  implicitly  takes  into 
account  bar  formation. 


(b)  The  change  of  volume  due  to  the  perturbation  and  departure 
from  the  equilibrium  profile.  This  departure  is  characterized  by 
a  change  of  slope.  Therefore,  the  corresponding  variation  of  sand 
volume  (area  ACC  in  Fig.  3)  is 


^2  i  «»S  y 

dt  2  dt 


3S 


y 


at 


in  accordance  with  the  previous  assumption.  Therefore,  the  total 
variation  of  sand  volume,  dV/dt,  is  (adding  eqs.  4  and  5): 


dV 

dt 


(B  .  U.)  ^2.  .  (ya  -  y},)  f 


q  at 


(5) 


(6) 


This  variation  of  volume  is  due  to  the  variation  of  littoral  drift  along 
the  ox-axis  and  the  onshore-offshore  motion.  The  following  terms  are  included: 


(a)  The  discharge  of  sand  leaving  the  beach  per  unit  of  width 
includes : 

(1)  Qyy  due  to  loss  of  sand  by  wind. 

(2)  Qy0(x)  due  to  the  quantity  of  silt  contained  in  the 

bluff  and  which  tends  to  move  offshore  by  suspension.  This 
loss  occurs  only  in  case  of  erosion  (Dy^/at  s  0)  and  is  equal 
to  =  KgB  3y,,/3t,  where  is  the  percentage  of  silt  in 

the  bluff. 


(3)  Q yj-  due  to  the  loss  of  sand  from  the  beach  by  density 
current  during  a  storm.  Qyr  is  a  function  of  the  size  dis¬ 
tribution  and  density  of  material.  A  beach  of  fine  material 
(<0. 1  mm)  tends  to  erode  more  rapidly  than  a  beach  of  coarse 
material  (> 1  mm) .  The  coarse  material  tends  to  move  along 
the  shore  while  the  fine  sand  moves  offshore. 
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The  determination  of  these  three  quantities  is  given  from  sand 
budget  investigations. 

(b)  A  general  term,  M(x,t),  expressing  the  local  variation  in 
the  sand  budget  due  to  loss  of  sand  by  rip  currents  along  groins, 
and  to  sudden  dumping  of  sand  in  case  of  beach  nourishment  or  flood. 

(c)  The  variation  of  littoral  drift  along  the  ox-axis  which  is 


°^0)  -  9a;(x  +  dx)  =  "  jT-  dx  (7) 

Many  formulas  in  literature  sources  express  the  rate  of  longshore  transport, 
Qy,  as  a  function  of  the  incident  wave  energy  along  a  straight  shoreline.  Long¬ 
shore  transport  is  also  a  function  of  the  sand  characteristics  (size  distribu¬ 
tion  and  density),  wave  steepness,  beach  slope,  etc. 

The  form  of  this  formulation  on  shoreline  evolution  is  of  paramount  impor¬ 
tance.  In  particular,  determination  of  the  relative  rate  of  sediment  transpor¬ 
ted  in  suspension  and  by  bedload  is  very  important  since  this  ratio  influences 
the  loss  of  sediment  by  rip  currents. 

Such  evaluation  is  beyond  the  state-of-the-art,  and  any  improvement  would 
require  a  major  effort  beyond  the  scope  of  the  present  investigation.  Any  im¬ 
provement  in  the  longshore  transport  rate  formula  could  eventually  be  introduced 
in  t lie  model  at  a  later  date.  Therefore ,  it  is  assumed  that  the  rate  of  sedi¬ 
ment  transport  is  independent  of  density  mid  size  and  depends  solely  on  the 


longshore  energy  flux, 


by  the  empirical  relationship 
%  =  7.5  x  1G3  I’s, 


where  Q,.  is  in  cubic  yards  per  year. 

l’f  is  in  foot-pounds  per  second  per  foot  of  shoreline  and  is  expressed  by 
the  relationship  (II. S.  Army,  Corps  of  Engineers,  Coastal  engineering  Research 
('enter,  1977) 


-  ££l  •: 


T(H0K.,)-'  sin  2a, 


where 


refraction  coefficient  from  deep  water  to  the  line  of  breaking 
incept i on 

wave  period 
deepwater  wave  height 

angle  of  the  deepwater  wave  with  the  shoreline 
angle  of  breaking  with  the  shoreline 


Note  that  this  equation  expresses  implicitly  the  rate  of  littoral  drift 
along  a  straight  shoreline  in  terms  of  breaking  wave  characteristics  (except 
for  the  deepwater  wave  height,  II,).  I'he  equation  is  assumed  to  hold  in  case 


of  gentle  beach  curvature.  The  refraction  coefficient,  K^,  and  angle,  a^, 
can  be  determined  as  functions  of  the  deepwater  wave  characteristics  HQ,  T,  aQ 
(or  a b)  and  the  angle  of  the  shoreline  at  breaking,  3ys/3x. 


In  the  case  of  a  groin  perpendicular  to  the  shore,  at  x  +  -  «,  the  deep¬ 
water  wave  angle,  aQ,  with  bottom  contours  is  equal  to  the  angle  a  of  the 
wave  with  the  ox-axis  since  the  shoreline  has  the  same  direction  as  the  ox-axis; 
at  x  =  o  (i.e.,  at  the  groin),  the  shoreline  becomes  parallel  to  the  incident 
wave  crest  very  rapidly.  Therefore,  ctQ  =  o  and 


-tan" 1 


3ys 

3x 


x 


o 


(10) 


In  tne  general  case,  for  any  value  of  x 


a  +  tan 


(ID 


The  breaking  wave  characteristics  of  wave  height,  H^,  water  depth  d^, 
and  the  angle  of  breaking  a ^  can  be  obtained  from  the  deepwater  wave  charac¬ 
teristics,  HQ,  T,  and  aQ.  aQ  is  given  by  equation  (11)  in  terms  of  a  and 
3ys/3x  which  takes  into  account  the  curvature  of  the  shoreline.  The  following 
equations  permit  a  determination  of  H^,,  d^,  and  a^,  provided  the  bottom  contours 
are  parallel  along  a  wave  ray  (Le  Mehaute,  1961),  but  could  be  curved  along  the 
shoreline  (Fig.  5): 


y 


Figure  5.  Fffects  of  wave  refractions 
on  a  curved  beach. 
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On  the  other  hand,  since  (see  eq.  lb) 


3x  3x  \°-"5  +  5,5  L0) 


and  by  taking  into  account  equation  (11) 


(VJ: 


Finally,  by  inserting  equation  (20)  into  equation  (19) 

0.25  *  S.5  "%  ,2y<; 

“  ‘  1  *  (3yVix)2  3X2 

In  case  of  wave  diffraction,  the  wave  height  varies  significantly  along  a  wave 
crest.  Then,  the  previous  refraction  coefficient,  k^,  lias  to  be  replaced  by 
a  combined  diffraction-refraction  coefficient,  such  as  knK-,,  (A  diffraction 
current  in  opposite  direction  to  a  longshore  current  takes  place  at  a  distance 
from  the  end  of  the  groin.)  The  variation  3Kp/3x  is  much  larger  than  3K.--/3x. 
Therefore,  in  analogy  with  equation  (18) 


cos  2a, 


In  a  diffraction  zone,  is  due  to  the  sum  of  variation  of  shoreline  direc¬ 

tion,  tan'1  3yJ./3x,  and  because  of  diffraction  the  rotation  of  the  wave  crest 
around  the  end  of  the  groin,  5  (Fig.  6).  0  is  the  angle  which  has  the  end 

of  the  groin  as  apex  ;uid  extends  from  the  limit  of  the  shaded  area  to  the  con¬ 
sidered  location  defined  by  x.  Figure  b  shows  that  0  =  aL;  -  O'  and  0'  = 
tan-1  x/t,  where  4  is  the  length  of  the  groin.  Therefore, 


.-i  3y<j 


30  _  30' 
3x  3x 


tan'1  - 


(XAY 


and,  differentiation  equation  (23)  and  inserting  equation  (24) 


Cy%) 


_  3?'yg  + 1 
2  3x2  8. 


i  *  (VJ 


where  3K^/3x  is  the  coefficient  of  variation  of  wave  height  uue  to  combined 
effects  of  wave  diffraction  and  wave  refraction.  It  also  includes  the  effect 
of  diffraction  current. 

An  empirical  formulation  for  determining  the  combined  effect  of  diffraction 
and  refraction  is  more  suitable  to  quantitative  analysis  of  a  real  sea  spectrum 


Figure  fa.  Diffraction  zone  notation. 


than  more  exact  theories  of  wave  diffraction  which  are  valid  for  periodic  waves 
over  a  horizontal  bottom  and  are  represented  by  a  Fresnel  integral.  For  this, 
it  will  be  assumed  that  the  energy  travels  laterally  along  a  wave  crest  as  well 
as  along  a  wave  ray.  The  lateral  speed  of  propagation  is  assumed  to  be  equal 
to  the  group  velocity  of  a  periodic  wave  of  average  period.  (Actually,  since 
the  problem  is  confined  to  very  shallow-water  waves,  and  the  longest  waves  of 
the  spectrum  diffract  most,  the  limit  of  the  diffraction  zone  is  defined  by  an 
angle  such  as  the  velocity  of  propagation  of  wave  energy  along  the  crest  and 
is  simply  /gD,  where  I>  is  the  water  depth.)  This  lateral  transmission  of 
energy  results  in  a  decrease  of  wave  energy  from  the  exposed  area  to  the  shaded 
area,  such  that  (Fig.  7) 


(-<’) 


where  s  is  the  distance  along  a  crest  from  the  end  of  the  groin.  Figure  7 
shows  from  simple  geometrical  relationships  that 


therefo  re , 


IP 


cos  (a  .  ♦  ir>°)  J 
oA  =  -v  tan  (  !Ti°  -  t  ) 
of  =  ■  tan  | 45°  *  t  j 

..  /T 

cos  (.»  ♦  45*7 


'  tan  (45° 
j  ,  ....  >  -  ° 


<_  ) 

li*1  dx 

■  1 


(-") 


(JX) 
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B 
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In  the  case  where  the  long  groin  is  in  the  previously  defined  wave  diffrac¬ 
tion  zone  as  in  Figure  6,  it  is  assumed  that  the  wave  energy  which  reaches  the 
groin  is  absorbed  by  friction.  It  is  also  assumed  that  the  combined  effects  of 
diffraction  and  refraction  of  a  wave  spectrum  can  simply  be  represented  by  a 
sinusoidal  variation  of  wave  height  along  the  breaking  line  (Mobarek  and  Wiegel, 
1966)  (Fig.  7). 


Assuming  that  KD(x)  is  of  the  form  sin  k(x  -  xQ )  as  shown  in  Figure  7, 
such  as  at  A  =  o,  at  C  =  1,  then  to  satisfy  these  conditions 


—  (~ 
21  \t 


»°  +  a<?)) 


(tan  (45°  -  a0 )  tan  (45  ■ 

which  after  some  simple  trigonometric  transformations  can  still  be  written 


(29) 


k  =  Ji  cos  2ao 


(30) 


to  satisfy  the  equation  expressing  the  conservation  of  energy  flux,  such  as 
given  by  equation  (29),  then 


a2  = 


k2  m 


fl  cos  2«c 


it  cos  (45  +  a0) 


2  sin 


and  K^(x)  =  H(x)/H0  is  given  by 


//2  cos  (2a())  V /2  /tt  cos  (2<*  )  \ 

KTn(xJ  =  \2  ~ sTn'_r  ■  /  sin  \ - 41 -  tx  +  f  tan  <4S°  -  S)lj 


(31! 


I- 
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By  differentiation  equation  (31) 


/ /f  cos  (2a0)  \1/2  it  cos  (2a0) 
dx  \  2  sin  /  “~4l 


/ ir  cos  (2a  )  \ 

cos  ^ - — — _  [x  +  l  tan  (45°  -  o^)]j  (32) 

Inserting  this  value  in  the  littoral  drift  equation  (22)  permits  the  mathemat¬ 
ical  model  to  completed. 

In  summary,  the  variation  of  volume  of  sand  dV/dt  given  by  equation  (6) 
is  equal  to  the  sum  of: 

(a)  The  loss  of  sand  by  wind  Q^. 

(b)  The  loss  of  silt  Q./s  =  KdB  3yd/3t  where  Kg  is  the  percentage 
of  silt  in  the  bluff. 


(c)  The  loss  of  sand  by  density  current  during  storm,  Q .  and 

lo^s  by  rip  current  (or  brought  upon  by  nearby  river) . 

(d)  The  quantity  of  sand  dredged  or  (at  the  opposite)  deposited 
during  beach  nourishment. 


(e)  The  variation  of  littoral  drift  along  the  ox-axis  is  SQ^/Sx. 

In  a  refraction  zone  alone,  SQ^/Sx  is  given  by  equation  (18)  where  K,. 
is  given  by  equation  (17).  oy,  is  given  by  equation  (16)  as  a  function  of 
ao>  and  ao  bX  equation  (11),  with  respect  to  the  ox-axis.  Also,  3«r/3x 
is  given  by  equation  (21). 


In  a  diffraction  zone  30^/Sx  is  given  by  equation  (22)  instead  of  equa¬ 
tion  (18)  and  the  diffraction  coefficient  Kg;  by  equation  (31),  and  3K  /3x 
by  equation  (32)  for  a  wave  spectrum,  oq,  is  now  given  by  equation  (23)  and 
3a^/3x  by  equation  (25). 


Since  all  the  phenomenological  equations  have  been  established,  it  is  more 
convenient  to  express  them  in  dimensionless  form.  The  "general"  equation  ex¬ 
pressing  the  sand  budget  balance  can  still  be  written.  (The  loss  terms  have 
been  dropped  for  simplicity  and  can  easily  be  included  if  necessary.) 


( B  [  x ,  t  ] 


3ys 

-  ~  + 


(y  .  -  y.  ) 


dp  JAy  M) 

dt  +  2  n<-*  dt~  =  Ti 


(33) 


For  an  analysis,  this  equation  is  considered  in  dimensionless  form. 

f  =  length 
B  +  D , 


where  is  chosen  later 


t  = 


At 


(B,,  *  D.)3 


Let 

(34) 

(35) 


Q  =  K^.K?  sin  2  t. 


(36) 
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The  general  equation  is  thus  transformed  to 


B  +  D  3y„ 

_ _s  + 

b0  +  dc  at 


(.9 a  ~  hi 


dD  +  1  Do  dA? 
dt  2  B0  +  Da  dt 


a  sin  2at 

3x  2 


The  hats  will  be  dropped  from  all  variables  at  this  point  in  all  discus¬ 
sions  except  those  relating  to  observed  results.  Also,  the  subscript  s  for 
y  and  Q  will  be  omitted.  Hence,  the  general  equation  used  in  the  following 
is 


B  +  Da  3y 

B o  +  Dc?  3t 


-  ybi 


dD 

dt 


1  dAy 

2  Ba  +  Da  dt~ 


.  K2r2  sin  2a, 
a  U  .1  v 

6X  2 


Recall 


(37) 


cos  ar 


7?  cos  at 


K*  = 


where  a[;  =  a  +  tan  1  3y/3x  and  ci£,  =  function  of  a  ,,  as  f(a(.,)  (the  function 

•e- 

4  sin  2 a 7-  =  K?  cos  sin  a^, 

3  .  /  3f  .  .  \  3ao  .  3oiJ 

37  cos  “o  Sin  ab  =  ^cos  ao  cos  ab  sTJ  ’  sin  ao  sin  aZ>)  sT"  =  F(ao)  7T 

The  general  equation  (37)  then  becomes  (after  some  rearrangements) 


f  depends  only  on  V0/ \.Q  as  previously  shown).  Therefore, 

^  if2i’ 

Note 


B0  +  D, 


_3y  _ _ 

3t  B  +  D, 


F(“0)  - 

1  + 


(dv/  \ 

\  /  aXJ 


32y 

T-s-  +  R(x,y,t) 


-  ?x- 


( 38) 


where 


+  ^c?  8  ct  dD 

R(x,y,t)  =  B  +  F(a7  37  "  (>'o  "  >’bl  dt 


1  D^»  dAy 


3K 


D 


and 


2  B,  +  D5  dt  '  3x 

a  =  a(x)  in  the  diffraction  zone 


(30) 


The  above  equation  is  the  general  dimensionless  form  which  gives  the  time- 
dependent  sand  budget. 


The  general  equation  is  nonlinear  and  appears  to  be  impossible  to  solve  ana¬ 
lytically.  As  it  is  also  difficult  to  solve  computationally  in  the  most  univer¬ 
sal  case  (where  lake  level,  beach  slope,  and  bluff-berm  height  vary  as  functions 
of  x,  y,  t) ,  several  simple  cases  are  described  in  detail  which  indicate  the 
behavior  of  the  equations.  Although  numerical  results  are  presented,  detailed 
descriptions  of  the  numerical  methods  used  are  given  in  a  later  discussion. 
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The  simplest  examples  are  those  of  beaches  having  negligible  or  identically 
constant  bluff-berm,  uniform  beach  slope,  and  uniform  lake  level.  In  these 
instances,  the  general  equation  (38)  reduces  to 


—  =  F(o0)  - —L - -y  +  2Kd  —E  cos  aQ  sin  ab 

at  x  +  (Vax)  3x 


(40) 


subject  to  the  appropriate  boundary  conditions. 


Suppose  (as  shown  in  Fig.  8)  the  beach  is  initially  zero,  bounded  by  a 
breakwater  which  is  a  complete  littoral  barrier,  and  subjected  to  waves  having 
uniform  positive  deepwater  direction.  The  boundary  conditions  are 


—  =  -  tan  a  at  x  =  0 
3x 


ay 

—  =0  at  x  =  L 
3x 


where  L  is  a  distance  along  ox,  large  enough  so  as  not  to  be  influenced  by 
the  groin. 


The  shoreline  is  calculated  for  a  fixed  H0,  T,  at  intervals  of  equal  At 
(shown  in  Fig.  8)  for  selected  times.  The  general  shape  of  these  curves  is 
similar  to  the  type  formulations  obtained  by  Pelnard-Considere  (1956) . 


v 


Figure  8.  Initial  shoreline  and  incoming  wave  direction  (.A); 

shorelines  at  successive  times  of  l.'-.t  ,  2At ,  5At  , 
lOAt,  and  20At  (B)  (vertical  scale  exaggerated). 
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An  extension  of  the  above  example  is  shown  in  Figure  9  where  both  sides  of 
a  breakwater  are  considered.  The  initial  position  of  the  shoreline  and  wave 
conditions  in  deep  water  are  as  before  (Fig.  9).  The  boundary  conditions 
become 

=*  0  at  x  =  -L 
3x 


3y 

3x 

3y 

3x 


0  at  x  =  -0  (simplification  of  more  exact  condition) 
-tan  a  at  x  =  +0 


-^  =  0  at  x  =  +L 
3x 


The  uniform  depth  theory  of  Penny  and  Price  (1952)  is  used  as  an  approxi¬ 
mation  (not  substantiated)  for  diffraction  about  the  end  of  the  breakwater. 
The  shoreline  is  calculated  for  various  multiples  of  a  fixed  At  (Fig.  9,A). 


V 


Figure  9.  Shorelines  at  successive  times  of  5At,  lOAt,  and  15At  (A); 

shorelines  for  time  2At  (B)  (vertical  scale  exaggerated). 
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Allowing  the  wave  direction  to  alternate  between  +a  and  -uc  at  a  fixed 
rate  gives  an  interesting  variant  to  the  previous  example.  The  boundary  condi¬ 
tions  are  as  before  for  a  >  o.  The  results  are  shown  in  Figure  9(B)  for  the 
same  values  as  the  physical  parameters  used  in  the  problem  of  Figure  9(A).  Of 
interest  is  that  the  undulatory  patterns  of  the  shoreline  in  Figure  9(A)  dis¬ 
appear  in  Figure  9(B).  Hence,  diffraction-induced  undulations  in  a  natural 
shoreline  probably  rarely  appear  since  offshore  wave  climates  are  often 
multidirectional . 

The  numerical  scheme  used  to  generate  the  preceding  examples  was  based  on 
the  use  of  implicit  finite  differences.  Such  schemes,  whether  implicit  or  ex¬ 
plicit  (or  both),  are  commonly  used  to  efficiently  solve  parabolic  problems. 
However,  even  in  the  case  where  only  refraction  is  considered  (Fig.  8),  the 
boundary  condition 


~  =  -tan  a  at  x  =  0 
3x 

numerically  gives  a  solution  which  initially  may  not  conserve  mass,  i.e.,  the 
integrated  transport  equation 


may  not  be  satisfied. 

Unfortunately,  this  feature  is  unavoidable  for  most  such  schemes  (the  ex¬ 
ceptions  are  discussed  below)  as  the  following  demonstrates.  Figure  10(A) 
shows  an  initially  straight  shoreline.  In  any  finite-difference  scheme,  after 
one  time  increment  At,  the  shoreline  is  bounded  below  by  the  solid  shoreline 
in  Figure  10(B).  This  shoreline  lias  the  least  possible  area,  A,  where 

Ax2 

A  =  tan  a  (41) 

The  conservation  of  mass  equations  requires 

AtQ(L)  =  At  cos  a  sin  >  A 
Thus,  At  must  satisfy  the  inequality 


'  -  - -  v-'-j 

-  sin  ■•■h  cos-  a 

Since  the  accuracy  (and  in  explicit  schemes,  stability  as  well)  depends  on  the 
ratio  X  =  At/Ax2,  the  above  inequality  places  a  lower  bound  on  the  accuracy  of 
the  solution  which  may  be  unacceptable  in  practice.  The  finite-difference  form 
of  the  equation  for  the  conservation  of  mass  may  lie  incorporated  directly  into 
the  numerical  scheme.  In  this  case  a  solution  exists  which  is  similar  to  the 
previous  case  but  shows  a  small  erosion  throughout  the  reach.  For  engineering 
applications,  the  primary  quantity  of  interest  is  the  amount  of  s;uid  on  a  given 
shoreline.  Then  it  is  more  important  to  conserve  mass  than  to  satisfy  the  shore¬ 
line  boundary  condition  as  written  in  the  present  form.  The  general  equation  is 


,  » 
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Figure  10.  Initial  shoreline  and  incoming  wave  direction  (A);  computed  and 
minimal  shorelines  for  finite-difference  scheme  of  time  lAt  (B) • 


now  used  to  derive  an  equivalent  equation  for  the  transport  Q  which,  although 
subject  to  similar  numerical  problems,  will  satisfy  the  transport  boundary  con¬ 
ditions  exactly. 

In  a  situation  where  only  refraction  is  important,  the  general  equation 
then  becomes 


ay  3Q 

Dt  ”  3x 


where  Q  is  cos  a,  sin  a/,,  is  f(at,),  and  «t,  is  a 
tiating  by  x  gives 


tan  1  3y/3x.  Differen- 


3  3y  320 
3t  dx  3x2 


(43) 


The  transport  function  Q  can  be  considered  as  a  function  of  a,  which  may 
be  solved  for  or, ,  such  as  a  =  g(Q).  Thus,  the  transport  equation  (43) 
becomes 


t  lie  re  fore , 


but 


3 


—  tan  (a  .  -  a) 


3  32Q 

—  t  an  (g  (Q)  -  a)  =  — r 

^  V  1 


3gCQ) 

3t 


cos2 


32Q 

(g(Q)  -  «)  -4 

3  x*' 


3g(Q)  =  dg(Q)  n 

3t  dQ  3t 

30  cos2  (g(Q)  -  a)  32Q 


there  fore 


Assuming  a  solution  for  this  equation  is  known,  the  shoreline  y  can  be  calcu¬ 
lated  from  the  equation 


y(t,x) 


y(o,x) 


+ 


dt 


In  practice,  the  equation  for  Q  is  not  solved  in  the  above  form.  Implicit  in 
the  above  formulation  is  the  assumption  that  the  function  g  exists.  However, 
as  shown  in  Figure  11,  g  is  not  single-valued  if  the  maximum  range  of  the 
angle  aa  is  greater  than  approximately  41°.  This  difficulty  may  be  removed 
by  considering  the  equation  for  Q  and  y  as  a  system  subject  to  the  boundary 
conditions  for  Q.  Note  that 

cos2  (g(Q)  -  op  _  dQ_/r  +  /3y\2] 

dg(Q)/dq  dajl  \3x/J 

Hence,  the  equation  for  Q  becomes 

3Q  _  dQ  1  32Q 

at  daa  l  +  (3y/3x)2  ^  (45) 


Figure  11.  Transport  function  QCa^,)  =  cos  a0  sin  za0 
for  selected  values  of  z. 
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A 


This,  together  with  the  equation  3y/3t  =  3Q/3x,  is  solved  in  a  cyclic 
scheme.  One  possible  method  is  the  centered  Crank-Nicolson  type  implicit- 
explicit  scheme  discussed  below.  Suppose  y  is  given  for  all  x  at  a  given 
time,  t,  and  that  from  t  the  wave  climate  is  specified  by  the  (constant) 
triple  (a, 11^,1).  Let 


L(t,x) 


\ 


dQ _ 1 _ 

dao  1  +  (3y/9X)2 


At 


Ax2 


L(t,x)  =  an  approximation  to  L(t,x) 


Integrating  the  Q  equation  gives 


where 


At 

32Q 

32Q 

" 

Q(t  +  At ,  x)  =  Q(t,x)  +  — 

lu.x)  ■jjr 

+  L (t  +  At , x)  — - 

JX 

t 

t  +  At. 

32Q  _  Q(x  +  Ax)  -  2Q(x)  +  Q(x  -  Ax) 
Tx3"  Ax2 


(46) 


(47) 


Integrating  the  y  gives 


where 


y(t 


+ 


At 

At,  x)  =  y  ( t ,  x)  +  — 

30  Q(x  +  Ax)  -  Q(x  -  Ax) 
3x  2Ax 


3Q 

3Q 

3x 

3  x  j 

t  +  At 

(48) 


(49) 


These  equations  are  solved  numerically,  subject  to  the  appropriate  boundary 
conditions,  in  steps  using  the  cyclic  algorithm: 

(a)  Let  L(t  +  At,x)  =  L(t,x)  Vx. 

(b)  Calculate  Q(t  +  At,x)  Vx  subject  to  the  appropriate  boundary 
conditions. 


(c)  Calculate  y(t  +  At,x)  Vx;  calculate  L(t  +  At,x)  and  set  this 
equal  to  L(t  +  At,x);  then  calculate  new  Q. 

(d)  If  new  Q  compares  with  old  Q,  stop;  if  not  go  to  step  (c). 


Tests  with  this  scheme  have  shown  that  it  converges  to  its  limit  after  one 
application  of  step  (c). 

This  method  can  easily  be  modified  to  solve  t he  equation  where  both  dif¬ 
fraction  and  variations  in  lake  level  are  allowed;  i.e.,  m  is  the  slope 
(refer  to  eqs.  37  to  40). 


3y 

3t 


3  2 

— -  Kn  cos 
3x  D 


%  sin  ab 


2  db 
m  dt 


(50) 
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For  convenience,  let  Q  =  cos  u0  sin  a£ 


Q  =  k£Q 


(51) 


As  before,  referring  to  equations  (43)  and  (44) 

1  dg (Q)  3Q 


8  3y  3 

—  —  =  —  tan  (g(Q)  -  u)  = 
3t  3x  3t 

Also, 


1 


3u 


cos2  (g(Q)  -  a)  dQ 


3t 


cos^  (g(Q)  -  a)  3t 


'Q  ,  3Q 

—  =  K2  —  +  2Kn  — -  Q 

3t  D  3t  D  3t 

from  equation  (51). 

The  second  term  in  each  of  the  above  two  equations  is  negligible  in  physi¬ 
cal  situations  of  usual  interest  where  the  distance  between  the  shoreline  and 
the  tip  of  the  breakwater  is  large  compared  to  the  distance  the  shoreline 
changes  during  a  time  At. 


Therefore,  the  transport  equation  becomes 
3Q 


% 


at 


dg/n 


cos2  (g(Q)  -  a) 


32Q 

si2 


(52) 


and 


(53) 


3y  _  3^  _  _1  dD 
3t  3x  m  dt 

This  system  is  solved  using  the  same  type  of  algorithm  as  used  previously. 

In  the  present  situation  where  only  refraction  is  important,  several  approx¬ 
imations  are  possible  which  produce  problems  having  analytic  solutions.  The 
most  direct  approximation,  and  essentially  the  assumption  of  Penard-Considerc 
(1956),  is  to  approximate 


32y 


3y  1 

— —  =  (z  cos  an  cos  ah  -  s m  sin  a/)  - ; - ; — r?  — 

3t  i +  (3yAx)  3x 


(54) 


where  z  =  0.25  +  5.5  H0/L0  (see  eq.  16),  subject  to  the  boundary  conditions 


by 


3  y 

—  =  -tan  a 

3.x 

X  =  o 

y (x,  t )  =  0  for  x  = 


—  -  a  — 

at  ax2 
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i 


A 
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where  a  is  a  constant.  For  the  standard  breakwater  problem,  the  most  logical 
choice  for  this  constant  is  given  by 

2  1 

a  =  - — = —  =  (z  cos  a,,  cos  a;,  -  sin  a„  sin  a;,) - — - — rr  (55) 

1  +  tan2  a  1  +  (V^x)  Jx  =  o 

since  the  shoreline  in  this  case  is  principally  governed  by  its  behavior  at  the 
breakwater.  This  problem  (defined  by  eq.  55)  has  the  solution 

y(x,t)  =  2  tan  a  e"x  -  tan  a  x  erfc^£-*^  (56) 

which  is  the  same  as  that  of  Penard-Considere  except  that  the  constant  a  lias 
been  changed.  This  problem,  however,  doesn't  conserve  mass  since 

00 

- —  I  y(x,t)  dx  =  z  sin  a  cos  a  sin  z a  cos  a 
;>t  -() 

When  this  approximation  is  used  in  the  transport  equation  for  Q,  the 
problem  becomes 

3Q  _  32Q 

3t  3x^" 

subject  to  the  boundary  conditions 

Q(x  =  o)  =  0 


which  has  the  solution 


Q(x  =  °°)  =  cos  a  sin 


Q(x,t)  =  cos  a  sin  a:,  erf 


U.«) 


Integrating  the  equation, 

3y  = 

)t  Jx 


gives 


y(x,t)  =  cos  a  sin  <,  y  t1/2  e-x"/4at  _  A  Crfc f  * 

t  oia  '  \v4a 


which  is  of  the  same  form  as  the  previous  solution  (eq.  56). 

III.  INPUT  DATA 

1.  Description  ;tnd  History  of  Navigation  Project  at  Holland  harbor. 

The  input  data,  which  arc  pertinent  only  to  the  present  investigation,  deal 
exclusively  with  the  area  near  Holland  Harbor,  Michigan  (11. S.  Army  l-ng i nee r 
District,  Detroit,  1975).  Construction  at  Holland  Harbor  began  about  1856  when 
the  city  of  Holland  cut  through  a  narrow  tongue  of  land  between  hake  Michigan  and 
Lake  Macatawa  (Fig.  12).  The  present  dimensions  of  the  navigation  project  were 
established  in  1909.  'Hie  existing  navigation  structures  have  been  constructed. 
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Figure  12.  Layout  of  shoreline  reaches  characterized 
by  different  beach  conditions. 

reconstructed,  and  repaired  by  segments  during  different  periods  over  the  past 
117  years.  In  general,  the  north  side  shows  accretion  while  the  south  side 
shows  accretion  from  the  breakwater  to  about  1,200  feet  south  and  then  general 
erosion  farther  south. 

The  following  areas  (Fig.  12)  appear  to  be  affected  by  the  navigation  struc¬ 
tures  as  evidenced  by  aerial  photos,  condition  surveys,  and  plat  maps  for  the 
period  of  record  (1849  to  1945): 

(a)  North  of  Holland  Harbor  (Faglecrest  Road  to  the  North  Pier  of 
Holland  Harbor).  It  is  considered  that  the  accretion  fillet  north  of 
the  piers  has  been  relatively  stable  since  1933.  After  that  date,  the 
predominantly  southward-moving  littoral  drift  has  been  diverted  lake- 
ward  resulting  in  rising  nearshore  elevations  north  of  the  breakwater 
and  deposition  of  material  in  the  entrance  channel.  This  reach  of 
shoreline,  about  5,000  feet  in  length,  is  characterized  by  increasing 
accretion  from  north  to  south  (Fig.  12,  reach  1). 

(b)  South  of  Holland  Harbor. 

(1)  Holland  Harbor  entrance  channel  to  a  point  200  feet  south 
of  the  Ottawa-Al legan  County  line  (Fig.  12,  reach  2).  Reach  2  extends 
about  2,000  feet  south  of  Holland  Harbor,  and  consists  of  a  sand  beach 
about  50  feet  wide  backed  by  low  sand  dunes.  The  shoreline  adjacent  to 
the  south  pier  accreted  at  the  rate  of  9.6  feet  per  year  from  1871  to 
1944,  a  total  of  700  feet.  This  area  of  accretion  diminishes  progres¬ 
sively  southward  from  the  piers  to  a  point  200  feet  south  of  the  county 
line  where  the  shoreline  before  major  construction  (1871)  coincides 
with  that  of  the  1929  and  1944  surveys. 

(2)  From  200  feet  south  of  the  Ottawa-Al legan  County  line  to 
147th  Avenue  (Fig.  12,  reach  3).  This  reach  is  about  2,400  feet  in 
length  and  is  characterized  by  high,  stabilized  s;md  dunes  (up  to  120 
feet  above  low  water  datum)  which  are  being  undercut  by  wave  action. 
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(3)  From  147th  Avenue  to  a  point  approximately  1,900  feet  south 
of  146th  Avenue  (Fig-  12,  reach  4).  Reach  4,  about  4,500  feet  long, 
is  characterized  by  eroding  sand  dunes  up  to  225  feet  high.  Figure 
12  shows  that  erosion  is  greater  at  the  northern  end  of  the  reach  and 
decreases  progressively  southward.  At  the  southern  end  of  the  reach 
the  shoreline  before  construction  of  the  navigation  structures  coin¬ 
cides  closely  with  that  of  the  1929  and  1944  surveys.  Fro si  on  appears 
to  be  greatest  at  a  point  about  800  feet  south  of  147th  Avenue  where 
the  shoreline  receded  about  220  feet  from  1866  to  1945,  or  an  average 
of  5  feet  per  year. 

The  rates  of  movement  of  the  shoreline  near  the  harbor  have  decreased 
(1945-73)  but  the  trends  of  accretion  on  the  north  side  and  erosion  on  the 
south  side  have  continued.  The  present  lower  rates  of  shoreline  movement 
are  due  to  several  factors: 

(a)  The  rate  of  shoreline  movement  due  to  the  navigation  structures 
decreases  with  time. 

(b)  Local  interests  have  built  more  shore  protection  structures  as 
their  property  was  threatened.  The  consequence  of  these  structures  is 
to  reduce  the  apparent  rate  of  erosion  locally,  but  this  is  accomplished 
at  tiie  expense  of  steepening  of  the  offshore  bottom  profiles  and  in¬ 
creased  erosion  downdrift.  The  shoreline  protective  structures  now 
extend  almost  continuously  from  1,000  feet  to  about  4,600  feet  south 

of  the  harbor. 

(c)  The  bluffs,  currently  under  erosive  wave  attack  between  5,000 
and  7,000  feet  south  of  the  harbor,  are  very  high  so  that  erosion  rates 
measured  in  feet  per  year  are  low  but  the  corresponding  volume  rates 
(in  cubic  yards  per  foot)  are  high. 

The  shorelines  apparently  have  become  comparatively  stabilized  since  about  1953. 
2 .  (leomo  rphologv . 

Figure  13  shows  the  nearshore  bathymetry  determined  by  a  survey  in  April 
1973.  Beach  profiles  from  this  survey  and  one  in  May  1945  are  summarized  in 
Figure  14.  Three  or  four  lines  of  longshore  bars  are  prominent  most  of  the 
time  and  occur  in  water  depths  to  approximately  50  feet.  The  first  bar,  in 
depths  of  1  to  4  feet,  can  change  rapidly  from  storm  to  storm  and  is  often 
broken  into  short  segments  from  200  to  1,000  feet  long  by  rip  channels.  Adja¬ 
cent  to  the  breakwaters,  the  rip  channels  tend  to  elongate  lakeward  as  far  as 
the  second  bar.  The  second  and  third  bars,  occurring  in  depths  of  o  to  14  feet, 
appear  to  be  affected  only  by  larger  storm  waves  and  are  relatively  continuous 
along  the  shore.  Aerial  photos  show  that  these  bars  tend  to  lose  their  promi¬ 
nence  as  they  approach  the  breakwaters  from  both  sides,  suggesting  that  pre¬ 
vailing  erosion  exists  across  the  bar  crests  near  the  breakwaters.  lie  fourth 
bar  occurs  at  depths  of  15  to  2d  feet  about  1,500  feet  offshore  and  is  charac¬ 
terized  by  crescentic  crests  displaying  wavelengths  between  3,000  and  5,000 
feet.  A  crescentic  fourth  bar  witli  an  average  wavelength  of  about  5,000  feet 
existed  north  of  Holland  Harbor,  but  none  existed  on  the  south  side. 

Although  the  general  bathymetry  follows  the  coastline,  significant  departure 
from  this  trend  occurs  in  a  region  offshore  of  the  Holland  Harbor  breakwaters 
(Tig.  15).  The  24-,  2”- ,  and  29- foot -dept h  contours  meander  sharply  off  the 
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14.  Beach  profile  summaries  for  Holland  Harbor 
District,  Detroit,  1975). 


( f rom  ll.S.  Army  l.ngineer 
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harbor  entrance,  delineating  a  ndgelike  formation  which  extends  obliquely 
lakeward  from  north  to  south.  Deeper  contours  to  the  southwest  of  this  forma¬ 
tion  suggest  that  it  may  extend  as  deep  as  50  feet.  The  ridge  has  a  maximum 
height  of  approximately  7  feet,  and  occurs  due  west  of  the  harbor  entrance. 
Between  the  ridge  and  the  harbor  entrance  is  a  steep  trough  about  34  feet  deep. 

A  lack  of  survey  data  before  the  harbor  construction  prevents  a  reliable  con¬ 
clusion  as  to  whether  or  not  this  underwater  ridge  was  caused  by  the  presence 
of  the  breakwaters.  However,  the  possibility  cannot  be  ruled  out  that  it  was 
formed  by  the  littoral  drift  arriving  predominantly  from  the  north  and  flushed 
lakeward  by  the  breakwater.  This  material,  unable  to  return  shoreward  due  to 
the  absence  of  significant  swell  activity  in  this  region,  could  have  accumulated 
to  form  such  a  formation  on  the  offshore  bed  over  a  long  period.  Implications 
of  this  interpretation  are  significant  because  it  may  mean  that  little,  if  any, 
bypassing  material  across  the  harbor  entrance  could  be  expected  to  reach  the 
southern  shore  of  Holland  Harbor. 

Many  of  the  dunes  reach  a  height  of  150  feet,  and  in  a  few  places  exceed 
200  feet  above  the  lake  level.  The  highest  dunes  are  confined  to  a  belt  about 

1  mile  wide  but  lower  ones  occur  for  a  width  of  several  miles  near  Holland. 

Dune  building  by  wind  action  is  still  active.  North  of  Holland  Harbor  some 
low  dunes  are  actively  migrating  inland  or  leeward  from  the  bluff  onto  rela¬ 
tively  level  upland.  In  this  area  the  dunes  are  40  to  120  feet  high.  At  nu¬ 
merous  places,  slumping  of  vegetation  is  evident  on  the  bluff  slopes.  South 

of  Holland  the  dunes  are  higher,  frequently  more  th;in  200  feet  above  lake  level. 
The  dunes  are  well  vegetated  except  for  occasional  blowouts  on  larger  dunes. 

At  many  places,  dunes  descend  directly  to  the  lake  or  to  a  narrow  strip  of 
beach  up  to  approximately  50  feet  in  width.  The  only  area  near  Holland  Harbor 
where  a  substantial  beach  exists  is  about  2,000  feet  immediately  north  and  about 
500  feet  immediately  south  of  the  breakwaters.  Maximum  beach  width  in  this  area 
is  about  500  feet  on  the  north  side  and  150  feet  on  the  south  side,  indicating 
the  predominance  of  littoral  drift  from  t lie  north.  Berm  development  is  gener¬ 
ally  poor  for  many  miles  of  shoreline  north  of  Holland.  Berms  develop  more  fre¬ 
quently  south  of  Holland,  growing  to  an  average  height  of  1  foot  in  the  swash 
zone.  These  berms  arc  usually  truncated  at  positions  of  a  recessed  shoreline 
where  a  rip  current  has  gouged  a  deep  channel  across  the  surf  zone. 

5.  Littoral  Materials. 

The  dominant  littoral  material  which  comprises  dunes,  beaches,  and  nearshore 
lake  bottom  is  glacial  sediment  belonging  to  a  fine  sand  category  (.less  than 
<>.42  millimeter  in  median  diameter).  Scattered  pebbles  are  found  both  on  the 
dimes  and  the  lake  bottom.  Sand-sized  material  is  dominantly  quartz  (SO  to  85 
percent)  with  up  to  12  to  15  percent  feldspar.  Heavy  materials  represent  only 

2  to  5  percent. 

a.  North  of  Holland  Harbor.  An  analysis  of  bluff,  foreshore  (.above  low 
water  datum)  and  bottom  samples  indicates  that  the  majority  of  tnese  sands 
belong  to  the  fine  sand  category.  (The  Unified  Soil  Classification  System 
defines  fine  sand  as  ranging  in  size  from  0.42  to  0.074  millimeter  (1.25  to 
3.75  phi),  and  medium  sand  as  ranging  from  0.42  to  2.0  millimeters  or  1.25  to 
-1.00  phi).  North  of  the  harbor,  the  average  bluff  and  foreshore  sand  consists 
of  88. 0  percent  fine  sand  and  11.7  percent  medium  sand;  the  average  lake  bottom 
sand  consists  of  05.0  percent  fine  sand  and  5.0  percent  medium  sand.  A  general 
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trend  exists  for  the  sand  size  to  become  progressively  smaller  from  foreshore 
to  offshore.  Average  median  diameters  are  0.54  millimeter  at  the  foreshore, 
0.28  millimeter  at  depths  of  5  and  10  feet,  0.25  millimeter  at  depths  of  15  and 
20  feet,  and  0.  IS  millimeter  at  depths  of  25  and  50  feet.  Pune  and  bluff  sand 
is  somewhat  filler  than  the  foreshore  sand,  but  coarser  than  the  lake  bottom 
sand.  There  is  little  systematic  variation  in  sand  size  with  distances  from 
tlie  breakwater.  Sorting  ranges  between  0.2S  and  0.40  phi,  indicating  very  good 
sorting;  sorting  varies  little  between  dune,  foreshore,  and  the  5-foot  depth 
and  becomes  poorer  toward  deeper  water,  reaching  approximately  0.40  phi  at 
depths  of  20  and  50  feet. 

b.  South  of  Holland  Harbor.  The  majority  of  the  sands  from  the  dune,  the 
bluff,  the  foresnore,  and  the  lake  bottom  on  the  south  side  of  Holland  belong 
to  the  fine  sand  category.  Generally,  only  minor  differences  exist  between 
average  sands  from  the  north  and  the  south  sides  of  Holland.  The  average  bluff 
and  foreshore  sand  consists  of  S5.1  percent  fine  sand  (finer  than  0.42  milli¬ 
meter)  and  14.9  percent  medium  sand.  The  average  lake  bottom  sand  consists  of 
94.7  percent  fine  sand  and  5.5  percent  medium  sand.  Median  diameters  decrease 
from  foreshore  toward  offshore,  about  0.25  millimeter  at  depths  of  5  and  10 
feet,  about  0.20  millimeter  at  depths  of  15  and  20  feet,  and  about  0.  1  <■>  milli¬ 
meter  at  depths  of  25  and  50  feet.  Dune  and  bluff  sand  is  finer  than  the  fore¬ 
shore  sand  but  coarser  th;in  the  lake  bottom  sand.  The  lake  bottom  sands  are 
distinctly  coarser  immediately  south  of  the  harbor  and  also  along  lines  11  and 
12  (I;ig.  15),  about  (>,000  feet  south  of  the  harbor.  Sorting  displays  a  wider 
scatter  than  on  the  north  side  of  Holland,  r;inging  between  0.5  and  0.5  phi. 

As  on  the  north  side,  sorting  becomes  poorer  from  foreshore  toward  offshore. 

The  principal  sources  of  littoral  material  are  the  beach  and  the  bluff. 
Generally  along  eastern  Lake  Michigan,  and  particularly  near  Holland  Harbor, 
the  supply  of  littoral  material  from  major  inland  water  runoff  is  negligible. 
For  this  reason,  the  long-term  geologic  trend  appears  to  be  for  sediment  from 
shore  erosion  to  fill  in  the  lake. 

c.  Shoaling  and  Dredging.  Shoaling  and  dredging  records  indicate  that  the 
annual  shoaling  in  the  entrance  channel  of  Holland  llariior  averaged  about  25,000 
cubic  yards  for  the  period  1965-70.  The  pattern  of  accretion  shows  some  yearly 
variation  and  indicates  that  material  encroaches  into  the  entrance  area  from 
both  north  and  south.  In  addition,  the  relative  amount  of  drift  from  north 
and  south  appears  to  vary  from  year  to  year.  Figure  16  shows  conditions  at 

t he  harbor  entrance  before  dredging  in  1975  and  1974. 

4.  Meteorology. 

The  dominant  factors  are  winds,  waves,  and  water  level  variations.  Tne 
winds  and  waves  are  directly  responsible  for  sediment  movements,  and  fluctua¬ 
tions  of  water  levels  separate  the  regimes  of  the  two  areas  affected,  i.c., 
the  backsiiore  above  the  waterline  and  tne  foreshore  below  the  waterline. 

a.  Kinds,  hinds  are  a  dominant  force  affecting  the  Holland  llariior  area; 
they  produce  a  number  of  major  effects  by  (a)  generating  a  force  for  waves, 

(b)  causing  lake  level  changes,  and  (c)  transporting  sand  across  the  beaches, 
particularly  the  finer  sediment  sizes.  hind  data  (available  from  ship  reports) 
over  the  southern  half  of  Lake  Michigan  (south  of  44°  N.)  show  that  more  than 


21,000  observations  were  taken  from  1963  to  1973.  These  data  have  been  sum¬ 
marized  for  a  12-month  period  (Fig.  17, A)  and  for  a  9-month  ice-free  period 
(Fig.  17, B).  A  3-month  ice  period  would  correspond  to  a  severe  winter. 

b.  Waves .  Sources  of  wave  data  for  Lake  Michigan  were  (a)  Summary  of 
Synoptic  Meteorological  Observations  (SSMO)  (National  Oceanic  and  Atmospheric 
Administration  (N'OAA) ,  1963-73);  (b)  Saville  (1953);  (cj  Cole  and  Hilfiker 
(1970);  and  (d)  Littoral  environment  Observation  (LliO)  program  data  from  the 
U.S.  Army  Coastal  engineering  Research  Center  (CFRC) . 

Figure  18  shows  deepwater  wave  roses  for  an  average  12 -month  period  and  for 
an  assumed  9-montb  ice-free  period,  respectively,  in  southern  Lake  Michigan. 
Comparison  of  these  wave  and  wind  roses  reveals  a  close  agreement  with  the  wind 
and  wave  statistics.  Dividing  wave  heights  into  two  groups  ( i . e . ,  by  heights 
lower  and  higher  than  4.1  feet),  the  low  waves  occur  most  frequently  from  the 
south,  while  the  high  waves  occur  predominantly  from  the  north. 

Figure  19  shows  exceedance  probabilities  of  deepwater  waves  based  on  SSMO 
and  hindcast  data  of  Saville  (1953)  and  Cole  and  Hilfiker  (1970).  Only  those 
waves  which  occur  during  a  9-montb  icc-free  period  and  are  traveling  shoreward 
arc  included  in  the  statistics.  The  three  statistics  generally  agree  quite 
well  for  wave  heights  up  to  about  8  feet.  For  wave  heights  larger  than  8  feet 
(which  occur  only  about  1  percent  of  the  time),  a  discrepancy  between  the  three 
is  evident.  The  SSMO  data  show  a  less  frequent  occurrence  of  higher  waves,  but 
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5.  Hydrology, 


a.  Water  Level  Variations.  figure  21  summarizes  the  observed  mean  monthly 
lake  levels  since  1920.  Records  are  available  back  to  ISnO. 

b.  Currents.  Currents  in  the  Holland  area  are  variable.  They  generally 
appear  to  take  a  northward  set  during  the  summer  ;ind  a  southward  set  during  the 
winter  in  response  to  the  prevailing  win's  although  reversals  are  frequent. 
There  also  appears  to  lie  evidence  of  a  rotational  current  system  in  Lake 
Michigan  which  is  not  very  well  understood.  In  the  immediate  nearshore  area 

it  is  likely  that  the  currents  are  wave  induced  and  will  propagate  alongshore 
in  the  direction  of  the  dominant  waves  at  the  time.  Inshore  of  the  first  liar, 
rip  currents  are  extremely  well  developed  and  can  add  a  significant  amount  of 
offshore  water  movement  to  the  general  longshore  current.  Intensity  of  the  rip 
currents  is  evident  from  a  number  of  rip  channels  which  are  gouged  1  to  2  feet 
below  the  adjacent  ground.  Longshore  currents  approaching  Holland  Harbor  are 
deflected  lakeward  along  its  breakwaters.  Aerial  photos  indicate  that  the  de¬ 
flected  current,  after  bypassing  the  harbor  entrance,  will  move  a  considerable 
distance  lakeward  before  turning  parallel  to  ;uid  toward  the  shore  off  the  down- 
drift  coast. 


c.  1 ce .  floating  ice  is  common  in  the  Holland  Harbor  area  during  tno 
winter  months.  Ice  accumulates  on  the  beach  by  wind  and  wave  action  and  re¬ 
portedly  often  extends  in  a  solid  mass  from  the  shoreline  to  more  than  1,000 
feet  offshore.  V.'ave  heights  are  substantially  decreased  within  an  ice  field 
and  Ice  along  the  shore  blocks  wave  action  to  some  extent.  However,  the  action 
of  ice  can  accelerate  erosion  processes;  e.g.  ,  ice  being  pushed  upon  the  beach 
will  loosen  consolidated  beach  and  bluff  material  and  some  of  the  sediment  be¬ 
comes  embedded  in  the  ice  at  the  shore  and  over  the  longshore  bars  to  be  carried 
elsewhere  if  the  ice  drifts  away  during  the  spring  thaw.  The  average  ice  season 
at  Holland  Harbor  extends  from  late  December  to  late  March. 

b.  Littoral  Drift  listimate  from  Wave  Statistics. 

Littoral  drift  has  been  computed  from  the  available  wave  statistics  dis¬ 
cussed  earlier.  The  governing  relationship  used  was  (iJ.S.  Army,  t.'orps  of 
engineers,  Coastal  engineering  Research  Center,  1977) 

')  =  0.  mm  1 33  I,. 

where  Q  is  the  potential  littoral  drift  in  cubic  yards  per  day,  and  1.  the 
longshore  momentum  flux  in  feet  per  pound  per  day  per  foot  of  beach.  The  re¬ 
lationship  between  i  and  <  is  determined  from  linear  theory  in  assuming 
that  the  bottom  contours  are  straight  and  parallel.  The  number  of  ua  es  per 
day  for  each  wave  height  period-direction  class  is  given  by 

8(>4  ■■  1U-’ 
n  - - j -  p 

where  p  is  the  probability  of  occurrence  of  that  wave  class.  Computations 
were  made  using  the  wave  statistics  fr  m  SSMO  data,  Saville  (l9;o).  and  Cole 
and  llilfiker  (1970).  The  duration  oi  computation  is  for  t be  9-month  ice- free 
period  from  late  March  to  late  December. 
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Figure  21.  Mean  monthlv  lake  leve 


All  of  the  littoral  drift  computations  show  a  large  gross  drift  of  300,000 
to  500,000  cubic  yards,  defined  as  the  sum  of  all  north  and  south  movement, 
regardless  of  direction,  whereas  the  net  drift  is  much  smaller.  Individual 
yearly  predictions  show  a  wide  variability  in  net  littoral  drift  quantities. 

The  littoral  drift  is  predicted  to  be  from  south  to  north  with  a  mean  value 
of  more  than  265,000  cubic  yards.  The  SSM0  data  and  Coles  and  llilfiker's  wave 
statistics  lead  to  the  prediction  of  61,000  and  76,000  cubic  yards,  respectively, 
from  north  to  south.  This  direction  is  consistent  with  the  observed  growth  of 
the  fillets  on  each  side  of  the  harbor,  the  nortli  side  showing  a  much  greater 
accretion. 


Figure  22  summarizes  the  monthly  fluctuations  of  littoral  drift  based  on 
the  SSMO  data.  The  figure  shows  that  the  monthly  drift  is  large  during  late 
fall  when  the  activity  of  extratropical  cyclones  intensifies  in  these  regions, 
and  small  during  the  summer  months,  especially  in  June,  when  wave  heights  are 
generally  smaller. 


Figure  22.  Monthly  variations  of  computed 
gross  and  net  littoral  drifts 
for  Holland,  Michigan. 


7.  Analysis  of  Aerial  Photos. 

A  sequence  of  cigiit  sets  of  aerial  photos  of  Holland  Harbor  taken  over  the 
interval  27  July  1950  to  5  October  1973  was  analyzed  to  assess  the  long-term 
evolutionary  development  of  the  shoreline.  A  regression  analysis  was  nor formed 
to  correct  the  shoreline  for  lake  level  variations  and  to  determine  the  long¬ 
term  cros i on-accret ion  rates  along  the  coast. 
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Table  1  presents  shoreline  positions  (with  ;md  without  lake  level  correc¬ 
tions)  as  a  function  of  distance  along  the  coast,  as  determined  from  the  aerial 
photos.  1'he  corrected  shoreline  positions  are  referenced  to  a  common  lake  level 
of  578.5  feet,  as  determined  by  a  regression  against  lake  level.  figure  23 
shows  the  1975  Holland  Harbor  shoreline.  1'he  origin  of  the  coastline  coordinate 
is  taken  at  the  harbor,  positive  in  the  north  direction.  The  shoreline  position 
was  measured  at  290-foot  increments  along  the  const . 

A  regression  of  shoreline  position  against  lake  level  and  time  provided  an 
estimate  of  the  long-term  evolutionary  trend  of  the  shoreline  at  each  station 
along  the  coast.  The  regression  study  revealed  general  erosion  extending  south 
of  the  harbor  and  accretion  immediately  north  of  the  harbor.  Over  an  8,110- 
foot  span  south  of  the  harbor,  the  average  beach  erosion  rate  was  ().7\r>  foot 
per  year,  liver  a  4,0bl)-foot  span  immediately  north  of  the  harbor,  the  average 
accretion  rate  was  1.05  feet  per  year.  A  span  of  10,585  feet,  starting  4,930 
feet  north  of  the  harbor,  had  an  average  erosion  rate  of  1.28  feet  per  year. 
These  evolutionary  trends  of  the  shoreline  include  natural  effects  as  well  as 
any  effects  which  can  be  attributed  to  the  harbor. 

figure  24  shows  the  evolutionary  smoothed  trend  (defined  as  the  average 
accretion  or  erosion  rate  between  1950  and  1973)  along  the  coast  near  Holland 
Harbor.  Lroston  dominates  south  of  the  harbor,  accretion  immediately  north, 
and  erosion  again  farther  north  of  the  harbor.  The  smoothed  curve  in  figure  20 
demonstrates  that  the  accretion  result  in;,  from  the  harbor  ext  onus  about  5  ,i'(h) 
foot  north  of  the  harbor.  These  estimates  of  the  evolutionary  trend  of  the 
shoreline  are  subject  to  large  annual  and  spatial  fluctuations. 

8.  Comparisons  of  Profiles  from  1945  and  1975  Surveys. 

The  surveys  of  May  1945  and  April  1973  (lake  levels  at  578.4  and  at  580.1 
feet.  International  Croat  hakes  Datum  (K.'hD),  respect  i  vely)  permit  a  comparison 
of  profile  changes  from  which  sediment  volume  changes  can  be  estimated.  Since 
the  actual  profile  locations  in  the  two  surveys  did  not  coincide  (in  general), 
the  1945  survey  was  interpolated  on  the  lines  of  the  1973  survey  to  allow  a 
comparison.  Table  2  summarizes  the  cross-sectional  area  changes  and  the  volume 
computations  in  two  parts  (above  and  below  low  water  datum).  The  last  column 
in  the  table  gives  the  erosion  or  accretion  rates  determined  from  the  aerial 
pilot  os . 


The  average  volumetric  accretion  rate  on  the  north  side  (from  ground  survey 
data)  is  2.13  cubic  yards  per  year  per  foot  of  shoreline  for  1,950  feet  (450  to 
2,400  feet  north  of  the  breakwater)  compared  with  the  average  shoreline  accre¬ 
tion  rate  of  1.92  feet  per  year  for  this  same  region  determined  from  figure  20. 
The  average  volumetric  erosion  rate  on  the  south  side  of  the  harbor  from  the 
survey  comparisons  is  1.11  cubic  yards  per  year  per  foot  of  beach,  whereas  t  lie 
values  from  figure  20  for  the  soutn  stretch  (470  to  5,550  foot  south)  give  0.82 
foot  per  year  .as  an  average  rate  of  shoreline  erosion.  The  agreement  on  the 
north  side  indicates  that  the  approximation  of  1  square  foot  of  beach  per  cubic 
yard  is  reasonably  valid,  whereas  on  the  south  side  t h i s  is  less  true.  I  he  dis¬ 
crepancy  oil  the  south  side  is  related  to  the  bluff  erosion  contribution.  The 
main  erosion  rate  below  the  low  water  datum  level  is  only  ().!(>  cubic  yard  per 
year  per  foot.  lienee,  the  bluffs  contribute  at  least  an  additional  1.25  cubic 
yards  per  year  per  foot,  which  exceeds  the  conditions  for  the  approx imat ion  of 
1  square  foot  of  beach  to  correspond  to  I  cubic  yard  oi  material. 


46 


r.il'K'  1.  History  o('  short'  I  ine  positions  at  liollaiul,  Michigan, 
Iron  1  *>:>()  to  IP",  with  anil  without  lake  level 
a'  erect  >  ••ns.  1 

uitiiii  |  Ihitimr  idtifvifil  ( roa  bit(  lin# 

po* I t low  j  I  ft ) 

'  1  *  j  ;  1  If.g  S»|-I  iKsslniM,  1  I  ’trie  lw»ti  |~  19  Oct .  196:  1  5  Sept  ■  lyt1  |  »  Ki»  19oS  ]~S  Pel .  1973 

1  lata  (ivrli  _ 


1 — 

v  t»«’ 

**  ■  lull  *u 

1113  3b 

nth.  n 

118*. 68 

1187.04 

1178.86 

1146.32 

1 

-  .  0  l  V 

♦  4.  54 

*44  <4 

y«;  3* 

987.  31 

998. '5 

9bb. OJ 

941.45 

921. IS 

-  M*v  JO 

4  1. 4C 

•  33.0.* 

851.40 

851.40 

884 . 14 

835.07 

835.02 

75 7.  39 

1  ' .  Ot> 

Vv  S3 

'36.  8 

’69. S3 

'66 . 7b 

785.90 

744.9' 

755.  16 

716.41 

-i»-  I'd 

*.0,41 

6’1.. *4 

68'.  ft’ 

68  '.  6  7 

705.68 

6').  79 

b6 7.  b 7 

655.04 

-  U  .  UO 

6  - 1 . :* 

t  U  3b 

6 : :  i ' 

blO. 3b 

654.97 

637.00 

677.09 

624. 33 

-  S  '4  t  0 

*  \i  it 

6  38.53 

61  3  99 

613.99 

6  37.00 

6  IS. 90 

677.17 

624.33 

59’  M 

tv  '  .  8  v 

60S  .  HO 

597.6) 

599.75 

605. 80 

589  43 

585. 39 

S’ 3. 05 

SO*. 3b 

S.'3.94 

513. '5 

540. 31 

57  3.94 

515. ’5 

511. 7* 

,0 

s;  s.  44 

4’4.8. 

491. 19 

48'  00 

499. 38 

49).  19 

4  74.82 

460.57 

- <*u.  Ok 

4  9  3 . 00 

4  33.  •* 

4.*5.  '0 

447.0’ 

447. O’ 

447.07 

435.52 

388.93 

•V9>  I'O 

450.2* 

39.'  - 95 

S'6.51 

473. '0 

409.  37 

447.0' 

477.42 

317.21 

4. 

-V»J*o  o0 

412.60 

3.'  V  46 

3: ’.46 

384. ?' 

343.83 

360.71 

397 . 95 

307.05 

„  | 

-  <  *  \  .  00 

Vb.59 

5i9. :* 

311.09 

3 '6. 38 

343.  83 

37 '.46 

343.83 

296.11 

-4  >b0  00 

3’b.58 

34  3  k 3 

335.65 

J7o. >8 

35  3.  6b 

330. 73 

357.07 

32  7  .  52 

-•<>  .S' 

i’t.M 

s-.o: 

343. 83 

3)6.58 

360.71 

336.93 

335.65 

286.51 

-«<■«.  t»0 

5'J.  50 

34  3.  83 

34  5 .  8  3 

3*6.40 

409.  37 

380.40 

348. 74 

766.11 

£  fe- 

-  44  So.  00 

ki.4* 

540 .  5b 

34  3, 8  3 

J*-7. 93 

401. 14 

381.49 

343.83 

276.34 

-  .  ou 

J’t. 58 

i’6. S8 

.368.  39 

3’6. 38 

37t. 38 

376.38 

357.02 

307.05 

-55  10.00 

384.  ” 

3bC .  .*  1 

36 8.  39 

384  .  " 

396. 73 

39-  .  95 

368  .  39 

276.34 

-sano.uo 

39.'  .95 

36  0.  :i 

343.83 

3'c. 58 

407. ’b 

409.  37 

376.58 

307.05 

-♦>0*10.00 

4i:.t>o 

34  3 .8  3 

376.58 

384.  ” 

409. 37 

409 . 37 

368. 39 

307.05 

-6  *S7.  OU 

4.'  ..  ’0 

34  3.  8  3 

.368.  39 

384.  ” 

475. ’0 

397.93 

J’6. 5i 

276.  34 

-M  '0.00 

419. is 

357.02 

343.83 

397.95 

474.87 

409.  37 

392.95 

286.51 

-64bt\0U 

599. SO 

360 . 2 3 

35  2.02 

389  68 

419.15 

41'. 51 

41'. 51 

276.34 

.  r.-i.o  mi 

184.68 

368. 39 

3 '6. 58 

397.95 

3  ’  3 .  30 

409  37 

4C9  37 

307,05 

'  -40.00 

l’b.58 

.36  8.  39 

357.0.' 

3b 8.  39 

360 . 7  1 

396 . 7  3 

3'6.58 

307.05 

•*HV  OJ 

343.  *3 

3o8 . 39 

37'  46 

'60. 7  1 

3’6. 58 

3'9  85 

345.85 

786  58 

•ll.'C.OO 

3;  • .  *t 

360. 2  1 

79'. 99 

•>7'.  40 

786.5 3 

34  3  8  3 

370  91 

745.64 

-»♦  10.00 

i:?.  46 

34  3  .  83 

794. '1 

311.09 

783.  75 

37’. 46 

794. 71 

255.17 

-•VO.  00 

s:’.4b 

1498. 13 
139y. 89 

311.09 

1555.43 

794. '1 
1449.01 

770.  15 
lb  17  '4 

794.71 
1540.  70 

7 '8.  34 
1547.  54 

245.64 

1498.40 

*’0.00 

1  34  .  59 

1309.84 

1391. '0 

13‘  8.96 

1474.4s 

1 358.9b 

1 3 ' S . 33 

1310.08 

■l 

1 160.00 

12  36. 1 6 

1318.03 

1 3uJ . 79 

1375. 33 

l. JO.  19 

1 3U4 . 9  3 

1228.20 

Ujo.00 

1  -  3B.  16 

1170.6' 

1719  '9 

1736. lb 

1309.84 

17  74.  V 

1700. ’7 

1136  08 

1  '40  00 

j:o3.4. 

1179. 74 

1187.04 

1195.73 

17(  0.  '7 

1706  33 

1719.  '9 

1166.79 

st 

si 

2050. 00 

1 1 78 . Bb 

1090.44 

1154.30 

1195.73 

1711.60 

1 19b. 8' 

1711  60 

1115.61 

2120.00 

114b.ll 

1096.99 

1105.  18 

1167  4  8 

1744.  33 

1193.73 

1195.23 

1045.97 

•Ml).  00 

1113.36 

1096.99 

1088.80 

113’. 97 

1777.9' 

1135.93 

1167  48 

1043.97 

400 .  Oil 

1074. O’ 

1064  .  .’4 

1090.44 

1105  18 

1146.  11 

114b.  11 

11 3’. 97 

1023  SO 

_ iiipjjio. 

u.M.-.u 

1015. 13 

_ KU.4  74 _ 

...lo?  3.89 

..  _ 

1 09b. 99 _ 

K  S-j  67 

Corrected  lata  level* 


0  } } 

*"  j  •  j  j  >  j 

T~ 

r~  /  j  /  ♦ 

1 1  '8  4 

j lb*  1  i 

1189  73 

-2  90 . 00 

974.29 

9*6, '2 

982.  38 

1 2  r£? 

986.43 

}  1  J  “ 

9 ' 3 . 46 

95'. 91 

931  33 

961.62 

-580.00 

8!'. 42 

84'.2b 

851.40 

8S8  38 

84  0 . 4  ’ 

*2  1  05 

81'. SS 

82  7  27 

-870.00 

761.86 

743.50 

769.5  3 

VO.  09 

761.93 

73’. 30 

*43  57 

’54.81 

-lltO.OO 

714. OS 

67b. 86 

6#;. e? 

690.84 

685.81 

664.94 

ft'9.72 

686  83 

-1450.00 

667.60 

633.59 

637.21 

643. J’ 

62  8 . 30 

622.47 

642  81 

-1740.00 

629.99 

634.*’ 

Mi.99 

614.  1? 

630.84 

6.4  S3 

671  V 

626  19 

-2030.00 

595.91 

607.29 

605.80 

508.47 

593.93 

604  1  0 

58’. 30 

391  91 

-2320.00 

SbS.’S 

511.33 

52  J.P4 

51'. 90 

526  85 

S  1 9  05 

510. 37 

533.28 

-7610.00 

sia.ti 

479  4  7 

491. 19 

485.66 

4*2. 'S 

48'. .  87 

461  17 

48’  18 

-2900.00 

4 ’5. 24 

440.68 

475. '0 

445  .  *>S 

41'  *0 

4  34.30 

475.81 

427. r 

-3190.00 

4  3o.  12 

4CS.33 

376.58 

4  32 .  7’ 

365.  13 

42’  93 

404  ’5 

387  99 

-34  80.00 

404. 30 

334.67 

32  7.46 

388  89 

318.0? 

.33  1 . 96 

387  65 

348  J7 

X 

-377Q. O0 

3'0 . 06 

324.98 

311.09 

379.84 

37.3.45 

320  94 

335  68 

329  4  2 

•>  « 

-4060.00 

373. 38 

546.63 

335. 6S 

3'8. 18 

34  3.6  7 

32  7.54 

34  8.02 

343.50 

-4350.00 

3o7. '0 

359  79 

54  3.8  3 

581.07 

3.37  4  5 

348  OS 

174  54 

111  00 

K  % 

-4640.00 

356. '8 

358.29 

343.83 

394.67 

357  68 

369  #8 

378  09 

348  74 

K-  •*- 

-4930.00 

34').45 

352.84 

34  3.83 

399  9? 

357.78 

36  7  46 

376  79 

346  51 

-5270. 00 

3o3.  52 

383.63 

368 . 39 

380  6  1 

351.40 

368. S 7 

34  1 . 93 

347  U 

-5510. OQ 

369.84 

373.26 

368.19 

392 . 2 3 

349  59 

3’8  03 

349  *4 

330  96 

-3800.00 

330.42 

371  17 

34  3.83 

382.84 

36  3  62 

396  #0 

360  92 

369  70 

-6090.00 

3'1'J.OS 

353.69 

.376.58 

391.54 

.366  98 

395  ’8 

131  46 

574.80 

-638C.00 

406. ’7 

360  40 

3b8.  39 

394 . 2 J 

366. 55 

3'4  0: 

352  92 

370  99 

-6670.00 

397.  14 

3' 1.28 

34  3  *3 

403  96 

406.03 

38'  31 

365  44 

396  64 

-6960  00 

330.76 

37'. 04 

352.02 

399  30 

359.07 

398  2  7 

39  3  46 

372  33 

-’2S0.00 

3 '5.97 

377 . 80 

376.58 

398  33 

3  39  69 

398  5' 

JVS  88 

360  83 

-7540. UO 

363  26 

375. t  7 

357.02 

377.55 

334 . 20 

38  *  91 

3bb  11 

348  66 

-7830.00 

331. '9 

377.18 

327.46 

365.23 

54  5  .  7  1 

369  87 

331  79 

336  7’ 

-8120.00 

32  J .  83 

366  UO 

297.99 

330. 7* 

765.6; 

33  ?  7  1 

317  63 

2  78  77 

-8410.00 

32  3.0' 

347.6? 

794.71 

313.28 

269  54 

'73  O' 

789  23 

2  7t  87 

870.00 

1331.98 

1409.1’ 

1391  ’1 

1 3b4 . 26 

1391  30 

l>4n  35 

!  36 2  0T 

1363  12 

1160, no 

12?C.*»0 

I 322  42 

1318.03 

1310.4# 

1330  39 

12 ’5  81 

1786  93 

1300  10 

5 

1450.00 

121’.  10 

1257.84 

1219. '9 

1248.69 

1750  :i 

12-5  64 

17  3*  90 

12  31. 18 

JS  ^ 

1740.00 

1 192  1 3 

11*0  02 

1187  04 

1700.57 

119-  64 

1206  4  3 

1220  21 

3  ■£ 

2030.  **0 

in.:.  *4 

1141.  13 

1154.30 

1201. 74 

1170.91 

1  18  3  8'- 

1193  33 

1)80  72 

s  ■ 

2  320.0,1 

1121  49 

1111.99 

1105  18 

1 1  '4  .  ’9 

1167  40 

ll'O  6* 

1164  45 

*  16  ?  08 

o  2 

26  10  00 

10-V  '5 

1115.03 

1088  *0 

114  8  .'3 

116  3  55 

113.  37 

1136  71 

114’  05 

2900.00 

lOM  9' 

1110.21 

1099  44 

1117.  V 

1098  9 1 

1131  Cl 

1119  04 

109902 

3190. UO 

1018. j 1 

IU’5.61 

10*  1.75 

1090  39 

10*5  *b 

10*4  00 

10M  38 

1067  99 

_ 34 80. OQ 

1020.92 

1071.68 

108f'  09 

104  8  59 

1C  16  01 

10"  6 3 

1043  63 

‘Determined  by  aerial  pftolo*. 


47 


able  2.  Summary  of  cross -sect  loual  changes  aiul  computed 
vo  1  iiji'ot  r  i  c  changes  for  1  ■_>  1 7»  and  l9-3  survov-..  ■ 
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As  demonstrated  previously  in  the  analysis  of  the  aerial  surveys,  the  har¬ 
bor  influence  on  the  north  side  extends  about  -1,2(10  feet.  farther  north  it  i- 
assuined  that  the  observed  erosion  rate  will  he  close  to  the  natural  erosion 
rate  for  this  area  of  the  hake  Michigan  coast.  At  1,200  to  la, 400  feet  to  the 
north  of  Holland  Harbor,  the  average  erosion  rate  during  the  1950-73  period 
(from  aerial  photos)  was  1,28  feet  per  year  per  foot  of  shoreline.  If  the 
approximation  ot  1  cubic  yard  of  volume  per  square  foot  of  beach  is  used,  the 
corresponding  sediment  volume  loss  is  1.28  cubic  yards  per  year  per  foot  of 
shoreline  plus  the  erosion  rate  of  the  dimes.  The  bluffs  which  have  an  average 
height  of  about  -It)  feet  in  this  area  would  contribute  an  additional  sediment 
loss  of  40  1.28/27  =  1.90  cubic  yards  per  year  per  foot,  lienee,  the  total 

sediment  loss  to  the  offshore  in  this  area  is  estimated  to  he  averaging  l.'.lu 
♦  1.28  =  3.18  cubic  yards  [icr  year  per  foot  of  shoreline  during  the  past  23 
years.  Th’.s  rate  of  loss  is  considered  to  represent  the  natural  rate  of  loss 
for  littoral  sediments  in  the  Holland  area  due  to  all  factors  except  t  lie  navi- 
ion  st  ructurc ,  It  must  he  remembered  that  the  erosion  rate  of  3.1s  cubic 

yards  per  year  per  foot  of  shoreline  is  an  average  rate  for  the  past  23  year 
(1930-73).  Since  lake  levels  have  been  much  higher  during  the  past  few  year--, 
the  present  erosion  rate  must  he  expected  to  he  higher  than  the  average  and 
therefore  will  probably  exceed  3.18  cubic  yards  per  year  per  foot  of  he.Kh. 


9 .  Sediment  Budget. 

a.  yorth  Side.  An  independent  evaluation  of  the  net  littoral  drift  ha 
been  obtained  from  various  surveys  as  suinmur  i  ceil  below  in  Wiiieh  a  sand  budget 
analysis  is  made  for  the  north  side  of  Holland  Harbor.  flu-  shoreline  changi 
shown  in  figure  23  provide  information  for  an  area  measurement  of  acc  is  t  i  on. 
An  approximation  for  volumetric  accretion  can  subsequently  be  made  iron  tin 
relationship,  wherein  1  square  foot  of  change  in  beach  surface  area  equal  ] 
cubic  yard  of  beach  material.  This  relationship  has  previously  been  shown  to 
be  reasonably  valid  for  the  north  side.  Subsequent  correct  ions  my  he  neves 
sary  because  of  the  dune-building  phenomenon. 
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The  area  changes  corresponding  to  the  waterline  shown  in  Figure  21  were 
determined  for  a  length  of  shoreline  extending  2,550  feet  north  of  the  north 
breakwater.  These  areas  are  summarized  in  Table  3.  A  correction  is  applied 
for  the  variation  in  lake  level  by  the  relationship 


Correction  in  1,000  square  feet 


45  x  lake  level  difference  *  2,550 

1,000 


which  follows  the  assumption  of  an  average  beach  slope  of  1  on  45  for  this  area 
as  determined  from  the  1045  surveys. 


Table  3.  Summary  of  shoreline  accretion  rates  north 
of  Holland  Harbor  (1840-1045). 


tear 

Measured  area  from 
breakwater  to 
2,550  feet  north 
(1.0 DO  ft2) 

Lake  level 
(S.Y.  datura) 

(ft) 

Lake  level 
di f fe rence 

(ft) 

Area 

correct  ion 

(1,000  ft2) 

Net  area 

(1,000  ft2) 

Rate 

(1,000  ft2/yr) 

1849 

714 

581.0 

.... 

.... 

714 

-  13.  3 

185b 

49S 

582.1 

♦  1.1 

♦  126 

621 

*51.9 

1871 

1,181 

582.9 

♦  1.9 

♦  218 

1,399 

♦  3.4 

1906 

1,462 

581.5 

♦  0.5 

♦5  7 

1,519 

♦13.7 

1953 

2,175 

578. S 

-2 .  S 

-287 

1,888 

♦13.4 

194S 

1,941 

580.  1 

•  0.9 

•  108 

2,049 

The  results  of  the  computations  (Table  3)  indicate  rates  of  accretion  vary¬ 
ing  from  3,400  to  51,000  cubic  yards  per  year  if  the  general  rule  of  1  square 
foot  of  beach  area  corresponds  to  1  cubic  yard  of  beach  material  is  applied. 

These  accretion  tes  should  be  interp reted  according  to  the  history  of  the 
construction.  (lie  first  breakwaters  were  constructed  in  185b-b0  and  trapped 
most,  if  not  all,  of  the  littoral  drift.  Since  the  accretion  of  51,000  square 
feet  per  year  during  the  period  1850-71  is  an  average,  then  the  corresponding 
littoral  drift  rate  must  have  exceeded  51,000  cubic  yards  per  year  (it  must  be 
assumed  that  by  1871  some  of  the  littoral  transport  was  passing  around  the  break¬ 
water  heads  into  the  navigation  channel). 

The  accretion  was  slowed  to  5,400  square  feet  per  year  from  1871  until  the 
harbor  breakwater  extensions  were  started  in  100b.  Therefore,  the  littoral 
transport  passing  around  the  north  breakwater  must  have  been  quite  high  from 
1871  to  190b;  also,  tlie  building  of  new  dimes  by  wind  action  probably  occurred 
during  tli is  period. 

From  100C)  to  1035  the  accretion  averaged  15,700  square  feet  per  year;  after 
1055  tiiis  rate  dropped  slightly  to  13,400  square  feet  per  year.  However,  some 
inland  dune  building  could  still  be  taking  place.  Measurements  from  the  1015 
survey  shew  an  accumulated  volume  of  157,400  cubic  yards  above  elevation  585 
(1015  datum)  on  the  north  side  of  Holland  Harbor.  This  volume  would  only 
account  for  about  1,000  cubic  yards  per  year  from  1871  to  1045  but  windblown 
sand  losses  farther  inland  would  increase  this. 

The  detailed  analysis  of  aerial  photos  taken  from  1050  to  1073  indicates 
an  average  accretion  rate  of  about  l.b5  feet  per  year  per  foot  of  shoreline  for 
the  first  4,000  feet  of  snoreline  north  of  the  harbor,  a  total  of  b,b()0  square 
feet  per  year.  Die  accretion  rate  based  on  a  comparison  of  the  1045  and  10~5 
hydrographic  surveys  for  profiles  1,  2,  and  5  in  Figure  10  was  about  1,800 
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cubic  yards  per  year  for  the  first  2,550  feet  north  of  the  breakwater.  The 
incoming  littoral  drift  supply  must  be  providing  the  observed  accretion  since 
there  are  expected  offshore  losses  of  3.  18  cubic  yards  per  year  per  foot  of 
shoreline  for  the  first  4,000  feet  of  shoreline  north  of  the  harbor,  and 
material  is  apparently  being  accumulated  in  the  dunes  in  Holland  State  Park. 

To  summarize  the  present  situation,  it  is  believed  that  about  01,000  cubic 
yards  of  littoral  material  is  being  transported  toward  Holland  Harbor  from  the 
north.  The  disposition  of  this  material  is  summarized  in  figure  20.  About 
12,700  cubic  yards  per  year  is  lost  offshore  (5.18  *  4,000),  0,600  cubic  yards 
per  year  is  accumulating  on  the  beach  face  on  the  north  side,  and  about  1,900 
cubic  yards  per  year  is  being  accumulated  in  the  dunes  of  Holland  State  Park. 
Although  the  exact  amount  is  unknown,  an  estimated  1,300  cubic  yards  per  year 
of  sand  is  lost  inland.  Hence,  it  is  calculated  that  about  38,500  cubic  yards 
per  year  of  sand  arrives  at  the  harbor  entrance  but  at  the  most  25,000  cubic 
yards  is  trapped.  It  is  likely  that  some  of  the  25,000  cubic  yards  dredged  in 
the  harbor  entrance  comes  from  the  south  side;  therefore,  at  least  13,500  cubic 
yards  per  year  is  lost  offshore.  Consequently,  it  is  concluded  that  the  harbor 
structures  cause  a  loss  of  materials  from  the  littoral  region  of  at  least 
13,500  cubic  yards  per  year.  In  addition,  the  25,000  cubic  yards  per  year 
that  is  dredged  is  dumped  offshore.  The  overall  result  is  that  the  south  side 
of  tlie  harbor  is  being  starved  of  01,000  cubic  yards  per  year. 


b.  South  Side.  figures  27  and  2S  provide  information  for  an  area  measure¬ 
ment  of  shoreline  accretion  and  erosion.  The  beach  area  changes  are  summarized 
in  Table  4  by  areas  north  and  south  of  a  section  located  1,200  feet  south  of  the 
south  breakwater  (an  arbitrarily  selected  station  which  appears  to  correspond  to 
a  stationary  shoreline  point  for  the  period  1850-1945).  The  table  summarizes 
the  corrections  in  lake  level  variation  about  an  average  offshore  beach  slope 
of  1  on  24.5  for  the  first  1,200  feet  south,  and  an  average  beach  slope  of  1 


Table  4.  Summary  of  beach  area  changes  south  of  Holland  Harbor 


on  12.1  for  a  length  from  1,200  feet  to  3,600  feet  south  of  the  breakwater. 

The  values  shown  in  Table  4  can  be  added  for  a  representative  length  of  shore¬ 
line  of  3,600  feet  to  yield  the  result  that  an  accretion  rate  as  high  as  48,500 
square  feet  per  year  in  the  period  1856  to  1871  decreased  to  a  rate  of  about 
6,000  square  feet  per  year  from  1871  to  1006,  then  increased  again  to  20,500 
square  feet  per  year  following  construction  of  the  breakwater  extensions  after 
1006,  and  finally  showed  an  average  erosion  rate  of  8,000  square  feet  per  year 
from  1933  to  1945.  This  last  value  corresponds  to  an  erosion  rate  of  0.75  foot 
per  year  as  determined  from  analysis  of  the  aerial  photos.  However,  it  should 
be  emphasized  that  the  recently  observed  erosion  rates  reflect  two  effects: 

(a)  The  length  of  shoreline  near  the  harbor  has  been  protected  by  rubble  and 
groins,  and  (b)  the  unprotected  shoreline  farther  south  is  backed  by  dunes  witli 
an  average  height  of  about  120  feet  with  peaks  more  than  200  feet  above  the 
lake  level. 

The  rate  of  accretion  at  times  on  the  south  side  may  seem  somewhat  surpris¬ 
ing,  but  it  must  lie  recalled  that  the  littoral  drift  often  reverses  (.see  Tig. 

22)  and  in  some  years  may  he  completely  reversed.  The  variable  direction  of 
littoral  drift  is  also  evidenced  by  the  two  surveys  in  figure  12,  which  shows 
the  harbor  entrance  shoal  from  the  south  side  as  being  apparently  larger  than 
that  from  the  north  in  1973. 

Using  the  approximation  of  1  cubic  yard  of  sediment  per  square  foot  of  beach 
the  1950-73  observed  average  erosion  rate  for  about  9,000  feet  south  of  Holland 
Harbor  is  0.75  cubic  yard  per  year  per  foot  of  shoreline,  plus  the  bluff  contri¬ 
bution  of  0.75  x  120/27  =  3.33  cubic  yards  per  year  per  foot  for  a  total  loss 
rate  for  unprotected  shoreline  of  4.08  cubic  yards  per  year  per  foot  of  shore¬ 
line.  When  compared  with  the  erosion  rates  observed  for  the  shoreline  about  1 
mile  or  more  north  of  Holland  Harbor  of  5.18  cubic  yards  per  year  per  foot,  the 
difference  of  0.9  cubic  yard  per  year  per  foot  of  22  percent  can  be  attributed 
to  the  navigation  structure.  The  effect  is  very  small  and  should  not  be  ex¬ 
pected  to  be  readily  apparent.  The  rate  of  loss  of  0.9  cubic  yard  per  year  per 
foot  of  bluff  material  would  correspond  to  a  shoreline  erosion  rate  of  less  than 
0.2  foot  per  year  for  a  bluff  height  of  120  feet. 
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The  actual  loss  of  land  to  the  south  of  Holland  (attributed  to  the  naviga¬ 
tion  structure)  is  relatively  small,  at  the  most  22  percent.  This  figure  has 
been  derived  as  a  23-year  average,  but  with  current  high  lake  levels  the  pres¬ 
ent  natural  erosion  rate  is  certainly  higher  than  the  average. 

The  evolution  of  the  shoreline  at  Holland  was  studied  using  the  present 
model.  The  relevant  physical  data  and  the  estimates  of  offshore  sediment  losses 
were  used  in  the  analysis.  Interpolation  of  values  shown  in  Table  1  gives  the 
shoreline  every  100  feet  along  the  base  line.  The  channel  at  the  harbor  en¬ 
trance  was  collapsed  to  zero  so  that  the  south  reach  ended  at  0-  and  the  north 

at  0+ .  These  shorelines  indicate  that  the  shore  is  stable  at  the  breakwater; 
therefore,  when  the  direction  of  transport  is  toward  the  breakwater  it  is  as¬ 
sumed  that  the  sand  transported  to  the  breakwater  is  entirely  lost  offshore. 
Monthly  lake  levels  were  taken  from  Figure  21.  The  results  of  these  computa¬ 
tions  are  not  given  since  they  are  essentially  previously  described  results. 

The  average  beach  slope  at  the  waterline  was  determined  to  be  1:10;  the  beach 

profiles  in  Figure  17  near  the  breakwater  show  that  this  is  a  reasonable  esti¬ 

mate  although  the  slopes  are  not  constant  from  profile  to  profile.  The  height 
of  the  berm  is  assumed  to  be  10  feet.  The  depth  to  no  sediment  motion  was  esti¬ 
mated  at  30  feet,  based  on  visual  consideration  of  the  offshore  bathymetry  and 
the  use  of  Weggel's  method  (J.R.  Keggcl,  personal  communication).  An  offshore 
loss  of  3.2  cubic  yards  per  year  per  foot  of  beach  is  also  included  (see  Fig.  20), 

The  transport  equation  for  this  situation  can  be  written  (see  eq.  55) 


IT  _  1__  rifZm  _  I  dD  _  dL 
3t  3x  ^  U*'  m  dt  dt 


where 


30  3  ? 

Q  =  cos  a0  sin  *b  and  (KpQ) 

m  =  beach  slope  at  waterline 
D  =  dimensionless  lake  level 

■jj-  =  dimensionless  offshore  line  loss 

This  equation  will  be  applied  on  each  side  of  the  breakwater.  When  the  in¬ 
coming  wave  direction  gives  drift  toward  the  breakwater,  the  boundary  condi¬ 
tions  expressed  in  Q  are 


-  cos  a  sin  a, 


Conversely,  when  the  incoming  wave  direction  gives  drift  away  from  the  break¬ 
water,  the  conditions  become 


=  COS  t  SMI  <: 
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Choice  of  wave  climate  is  the  remaining  input  parameter  to  be  determined,  and 
is  the  most  controversial.  The  wave  climate  most  desirable  for  the  study  of 
shoreline  evolution  is  a  time  series  giving  wave  height,  H,  period,  i,  and 
direction,  D.  Unfortunately,  this  is  seldom  available  and  lienee  monthly  sta¬ 
tistical  summaries  must  be  used,  such  as  those  by  the  SSMO  for  southern  Lake 
Michigan.  One  possible  use  of  the  summaries  is  to  construct  monthly  times  t 
for  each  possible  (11,T,D)  triple,  i.e.,  t(H,T,D),  and  then  calculate  the  evo¬ 
lution  of  the  shoreline  as  the  (H,T,D)  triples  are  chosen  in  some  determin¬ 
istic  order  or  at  random.  This  method  would  be  computationally  very  expensive 
and  is  not  used.  The  most  simple  approach  is  to  assume  that  these  are  but  two 
(H,0,T)  triples  representing  the  gross  transport  north  and  south,  each  occur¬ 
ring  for  some  length  of  time  per  month.  The  entire  shoreline  is  alternately 
calculated  for  an  incremental  time,  assuming  the  direction  of  the  incoming 
wave  is  positive,  then  negative.  The  period,  T,  used  is  taken  to  be  the 
average  T;  i.e., 

T  -PHU'O  t 

1  -  lip  (11,  T) 


where  the  p(ii,T)  are  the  (II, T)  probabilities  given  in  Table  19  of  the  SSMO 
(National  Oceanic  and  Atmospheric  Administration,  1965-73).  The  choice  of  (II, 0} 
for  north  and  south,  denoted  (Ifyy, Dy)  and  (ll^.Dcj),  respectively,  must  now  be  made. 
This  choice  is  subject  to  the  condition  that  the  actual  northerly  transport,  as 
calculated  from  the  statistics  and  given  a  straight  shoreline  for  the  reach  of 
interest,  be  preserved;  i.e., 


t  .II-.’.T  cos  a  )  sin  t-,  =  t;,  ^  P ('*'(**)  PU<»U)  TIL  cos  x  .  sin  <. 

li,T  I)  giving 
north  transport 


where 


t;, 

P(T|U) 

p(!l,D) 

.i 

u, 

t  , 
II.. 

‘  ; 

I) 


number  of  hours  in  a  given  month 
conditional  probability  I  occurs  given  II 

probability  of  (11,1))  pair  (using  Table  IS  in  SSMO  as  a  data  base) 
I)  -  shoreline  orientation 

f  ('*,.) 

number  of  hours  the  average  wave  condition  exists 
average  wave  height 
I)  ■  -  shoreline  orientation 
a ve ra  ge  di rec f i on 


and  similarly  for  the  directions  giving  south  transport. 


The  average  directions  of  the  shoreline  at  the  breakwater  are  calculated 
from  the  historical  records.  The  directions  are  chosen  for  the  incoming  wave 
angles  since  the  complex  geometry  of  the  harbor  breakwater  slmlJs  the  nearby 
shoreline  from  waves  arriving  from  most  directions.  At  present,  the  time  dura¬ 
tion  of  waves  arriving  from  the  north  is  tttsuuicd  to  be  the  same  as  from  the 
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south ;  therefore,  t  =  t  •  =  0.f>  t^.  The  conservation  of  transport  equation  is 
then  used  to  calculate  the  average  wave  heights  lly  and  II  •  .  The  results  of 
these  calculations  are  given  in  Table  5. 


Table  5.  SSMO  averages 

for  Holland  Harbor 

Mont  h 

If. 

Hy 

T 

(ft) 

(ft) 

(s) 

Jan. 

1.3 

4.3 

6.  1 

lob. 

4.0 

3.  6 

6.0 

Ma  r . 

3.:> 

1.9 

6.3 

Apr. 

2.6 

2.  4 

5.5 

May 

2 . 3 

2.  1 

5.4 

Juno 

2.1 

l.S 

5.  3 

July 

1.8 

2.  1 

5.  3 

Aug. 

2.S 

2.5 

5.6 

Sept  . 

3.0 

2.6 

S.  8 

Oct  . 

3.  6 

3.  6 

5.9 

\'ov. 

1.0 

3.4 

6.0 

Dec. 

4.2 

3.8 

5 . 9 

Annua  1 

3 .  3 

2 . 9 

~  5.8 

•J*o  i*  L..  -  -SS;  L’.  -  *j7u. 


The  computer  program  described  in  the  Appendix  was  used  to  calculate  the 
evolution  of  the  shoreline  from  September  19<>7  to  May  19(>8.  The  historical 
1907  and  19o8  shorelines,  as  well  as  the  computed  1 9<>  <S  slice  lino,  are  shown 
in  figure  29.  The  calculation  assumes  that  Ax  =  100  feet  with  \  1  which 

gives  a  value  tor  At  which  varies  from  8  to  20  hours  depending  on  the  monthly 
wave  characteristics.  The  principal  discrepancy  between  the  predicted  and 
actual  1908  shoreline  occurs  near  the  breakwater.  Although  the  shapes  agree, 
there  is  an  erosion  in  the  calculated  shoreline  which  is  probably  due  to  the 
approximations  used  in  the  calculation  of  the  diffraction  coefficients,  and 
incoming  wave  angles  which  are  functions  of  x  in  the  shadow  region  of  the 
diffraction  zone.  The  unaltered  theory  of  I’enny  and  Price  (1911)  was  incor¬ 
porated  into  the  numerical  scheme  since  most  breakwaters  can  be  represented 
as  line  barriers,  and  hence  is  almost  always  useful.  However,  for  the  case 
of  Holland  Harbor  a  universally  valid  prediction  of  the  shoreline  would  re¬ 
quire  the  detailed  calculation  of  the  diffraction  effects  due  to  the  geometry 
of  the  breakwaters.  Also,  the  convenient  choice  o t  incoming  wave  direction 
obscures  the  fundamental  problem  of  how  to  properly  use  the  statistical  wave 
s  ununa  r  i  e  s . 


IV.  Cl i.MCl.IJS  1 0\S  AM)  HI. COM.' II  N DAT  IONS 

The  basic  idea  of  Pe  lna  rd-Cons  i  de  re  (19f>(>),  i.e.,  to  investigate  Mioreline 
evolution  by  concentrating  on  conservation  of  mass  as  a  special  one-dimensional 
problem,  has  boon  generalized  to  essentially  its  limits  of  app 1 i cab i 1  i  i v.  Ihese 
physical  processes  of  refraction  and  diffraction  (where  applicable)  have  been 
incorporated,  including  deterministic  variations  in  lake  level,  bluff  height, 
and  beach  slope.  The  inclusion  of  refraction  makes  possible  the  proper  use  of 
the  known  physical  relationships  between  wave  energy  and  littoral  drift  on  a 
priori  basis  without  necessarily  determining  these  as  results  from  the  past 
recorded  shorelines  at  a  given  locution.  Accurate  determinut ion  of  shore  be¬ 
havior  in  the  lee  of  a  breakwater  requires  inclusion  of  diffraction  in  some 
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*10(1  F  T 

Figure  29.  Computed  data  versus  actual  data. 

form.  This  could  be  done  in  a  heuristic  manner,  either  in  the  global  approxi¬ 
mation  described  previously  or  in  the  use  of  the  constant  depth  theory  of  Penny 
and  Price  (1952).  It  could  also  be  done  in  a  more  rigorous  manner  which  would 
include  the  effects  of  a  sloping  beach.  Thus,  quantitative  predictions  of  the 
shoreline  can,  in  theory,  be  attempted  in  situations  where  onshore-offshore 
transport  of  sand  is  either  negligible  or  is  known  from  other  sources. 

The  resulting  theory  is  presented  in  two  equivalent  forms,  one  in  terms  of 
the  behavior  of  the  shoreline  y(x,t)  alone,  the  other  expressed  explicitly  in 
the  longshore  transport  Q(x,t)  and  implicitly  in  y(x,t).  The  former  has  the 
advantage  that  numerical  schemes,  such  as  that  of  Crank-Nicolson  can  qualita¬ 
tively  indicate  the  behavior  of  the  shoreline  in  regions  of  rapid  change.  How¬ 
ever,  the  conservation  of  mass  is  difficult,  if  not  initially  impossible  to 
achieve  since  any  approximation  of  a  transport-derived  term  (i.e.,  a  term 
arising  from  0Q/)x)  will  alter  the  transport  balance.  The  later  form  allows 
the  use  of  analytical  or  numerical  approximations  in  the  transport  equation 
which  will  not  disturb  the  total  sand  content  of  the  system,  but  only  its 
local  distribution. 

The  most  severe  and  unavoidable  limitations  to  the  engineering  application 
of  these  methods  are  the  use  of  the  statistical  wave  summaries.  One  possible 
use  of  these  statistics  was  shown;  however,  many  others  are  possible,  lifficient 
and  accurate  use  of  the  offshore  wave  statistics  is  endemic  to  the  problem  of 
large-scale  shoreline  prediction,  and  must  be  achieved  before  any  theory  (whether 
one  line,  multiple  lines,  or  gridj  can  successfully  produce  accurate  results. 

The  problem  of  shoreline  evolution  sensitivity  to  time  step  in  the  input 
wave  climatology  would  require  further  research.  Despite  this  limitation,  if 
the  effects  of  wave  refraction,  wave  diffraction,  and  change  of  lake  level  are 
taken  into  account  (as  in  this  report),  and  if  the  method  is  generalized,  then 
a  mathematical  model  with  multiple  bottom  contour  lines  could  be  formulated 
which  would  (if  the  problem  of  wave  statistics  input  is  solved)  permit  a  cal¬ 
culation  of  the  evolution  of  tile  complete  bottom  topography. 
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It  is  important  to  point  out  that  wave  refraction  effect  on  shoreline  evo¬ 
lution  has  been  found  particularly  important.  It  is  particularly  necessary  in  J 

order  to  determine  a  planform  stability  criteria  which  can  be  established  from 
the  present  formulation. 

The  problem  of  shoreline  stability  needs  to  be  investigated,  both  physi¬ 
cally  and  numerically,  as  for  some  deepwater  wave  angle,  shoreline  perturbance 
may  increase  by  instability  instead  of  being  flattened  out. 

Perhaps  one  of  the  most  significant  results  of  this  mathematical  approach 
to  shoreline  evolution  is  to  point  out  the  need  for  more  research  to  quantify 
the  phenomenology  relevant  to  shoreline  evolution.  The  insertion  of  empirical 
parameters  has  allowed  the  investigator  to  fit  (to  a  large  extent)  a  form  of 

shoreline  evolution;  however,  the  processes  involved  in  each  of  these  param-  >• 

eters  are  not  well  known,  ;ind,  what  fits  at  Holland,  may  not  necessarily  apply 
elsewhere.  for  example,  the  research  topics  which  will  improve  the  model  are  ;  - 

(a)  onshore-offshore  movement;  (b)  quantity  of  sand  (silt)  lost  by  rip  currents 
or  density  currents;  (c)  percentage  of  sand  in  suspension;  (d)  distribution  of  ^ 

sand  discharge  as  a  function  of  the  distance  from  shore;  and  (e)  more  iupor-  j 

tantly,  how  to  treat  the  wave  climatology  in  finite  time  intervals  to  obtain  [' 

an  equivalent  result.  The  last  topic  may  be  the  most  difficult  since  the  noise,  1 

(i.e.,  daily  variation  and  effect  of  storm)  may  exceed  the  signal  (i.e.,  the  1 

long-term  trend).  Shoreline  evolution  is  due  to  a  succession  of  extreme  events 
separated  by  long-period  time  effects  of  equal  importance.  Therefore,  the  j, 

treatment  of  long-term  evolution  from  an  average  wave  climatology  is  question¬ 
able,  since  the  complete  time  history  may  have  to  be  considered  on  a  daily  basis. 

This  topic  can  be  investigated  by  a  sensitivity  analysis  to  wave  definition; 
e.g.,  the  present  model  could  be  used,  as  well  as  a  research  guideline  for  fill¬ 
ing  many  knowledge  gaps  in  shoreline  processes. 
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APPliNDIX 


COMPUTER  PROGRAM 

This  appendix  presents  the  listing  and  brief  explanation  of  the  computer 
program  written  to  investigate  the  behavior  of  a  shoreline  which  contains  a 
complete  littoral  barrier  at  x  =  o.  The  numerical  scheme  is  based  on  the 
finite-difference  method  of  C ran k  -N  i  co  1  s  on ,  used  to  solve  the  cyclic  nonlinear 
transport  equations  in  Q  and  y.  Although  the  program  is  written  expressly 
for  Holland  Harbor,  Michigan,  hopefully,  enough  explanation  is  given  to  modify 
the  program,  if  necessary,  to  suit  a  particular  application. 

The  program  consists  of  a  main  program,  which  is  simply  a  calling  routine 
to  the  controlling  subroutine,  and  eight  subroutines.  The  interrelationship 
of  these  programs  is  shown  in  the  figure  below;  a  brief  discussion  of  each 
subroutine  follows. 
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l-igure.  Program  structure. 


MAI. iQ  Calls  subroutine  TVOLVIIQ 

liVOLVF.Q  The  controlling  subroutine  which  organises  the  numerical  method  in  a 
global  sense  and  calls  various  other  subroutines  as  required.  The 
first  call  is  to  IMPHTO,  which  reads  the  controlling  parameters,  his¬ 
torical  shorelines  ;md  lake  levels,  ami  statistical  wave  information. 
For  each  month  in  the  time  period  of  interest,  subroutine  CAl.CkDq  is 
called.  (It  is  assumed  that  during  a  month  the  diffraction  coeffi¬ 
cients  change  negligibly.)  For  alternating  equal  incremental  times 
the  angle  of  the  incoming  waves  is  changed  ;tml  subroutine  SOl.VliQ  is 
called.  This  is  repeated  until  a  month's  time  has  been  completed. 
Then  this  is  repeated  until  another  historical  shoreline  is  reached. 
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INPUTO  This  routine  reads  the  information  required  to  define  a  particular 
problem.  Below  are  the  cards  required  with  an  explanation  of  the 
variables  required.  Title  cards  are  user  information  devices  which 
are  printed  but  never  used.  The  card(s)  are  numbered  by  their  order 
as  defined  by  READ  statements. 

CARD  1:  title  card 

CARD  2:  LAMBDA  ~  \  =  At/Ax2 

B  =  bluff  height  (feet) 

DO  =  depth  to  no  sediment  motion  (feet) 

DX  =  Ax  (feet)  ; 

BliRM  =  berm  height  (feet) 

BWL  =  length  of  breakwater  (feet)  : 

BWD  =  depth  of  breakwater  tip  (feet) 

ShOPli  =  slope  of  shore  profile  at  waterline  I 

CARD  3:  XLINli  =  offshore  line  loss  (square  feet  per  hour) 

♦ 

CARD  4:  title  card 

CARD  5:  N1  =  number  of  grid  points  to  left  of  breakwater 

N2  =  total  number  of  grid  points  on  both  sides  of  the 
breakwater 

Nl)  =  number  of  historical  shorelines 

NY  =  number  of  years  lake  levels  are  to  read 

INYR  =  year  of  the  first  historical  shoreline  (must  be 

year  of  first  lake  level  as  well) 

CARD  6:  for  all  years  with  historical  shorelines  read 

year  (I);  I  =  1,ND  i 

month  (I) ;  I  =  1 ,ND 

hour  (I);  I  =  l,Nl)  where  hour  is  the  day  *  24 
CARD  7:  for  each  grid  point  1  read  historical  shorelines  (feet) 

CARD  8:  title  card 

CARD  9:  for  each  year  read  lake  level  for  each  month 
CARD  10:  title  card  yj 

CARD  11:  for  each  month  (including  annual  which  is  not  at  present 
used)  input  (in  the  notation  used  for  Holland)  II  ,  !!■■, 

D„,  D.,,  T,  T. 

U  ll 

The  input  data  above  are  printed  without  explanation  along  with 
several  assays  used  internal ly.  Note  that  in  this  routine  i  is 
redefined  slightly  so  that  a  month's  time  is  exactly  2n.'.t  . 
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SOLVEQ 


CALCKDQ 

Q REG  ION 
S REG  ION 
CS 
CETC 
Note  I. 


Note  II. 


This  subroutine  calculates  the  shoreline  for  one  time  step  assuming  a 
previous  shoreline  (and  shoreline  angle  =  tan-1  ( ty/ <x) )  and  diffrac¬ 
tion  coefficients  as  calculated  in  subroutine  CALCKDQ.  The  scheme  is 
essentially  as  is  described  in  Section  2.  The  matrix  inversion  is 
achieved  using  a  standard  algorithm  for  tridiagonal  matrices. 

This  routine  calculates  the  diffraction  coefficient  and  incoming  wave 
angle  -a  =  x(x)  in  the  diffraction  zone  using  the  theory  of  Penny  and 
Price  (1952)  for  both  incoming  wave  directions  (the  angle  convention 
is  shown  in  subroutine  INPUTQ). 

Calculates  diffraction  in  the  "Q"  region. 

Calculates  diffraction  in  the  "S"  region. 

Calculates  complex  valued  Fresnel  integral. 

Calculates  a  given  depth  and  period. 

In  subroutine  EVOLVEQ 
LI  =  1st  year  index 
1.2  =  last  year  index 

Calculation 

begins  at  date:  year  (LI),  month  (LI),  hour  (LI) 
ends  at  date:  year  (1.2+1),  month  (1.2+1),  hour  (L2+1) 

In  a  general  program 
1.1=1 

L2  =  ND-l  where  NO  is  defined  in  subroutine  INPUTQ  Nl)  must  be 
passed  to  subroutine  EVOLVEQ. 

Program  as  given  in  this  report  is  set  up  to  allow  complete  bypassing 
(lost  offshore)  when  waves  are  not  diffracted,  i.e.. 


To  change  program  for  no  bypassing,  four  statements  must  be  changed 
in  subroutine  SOLVHQ . 

2  statements  after  1001 

change  QOLD  (N1P1)  =  (1-ID)  *  QOI.D  (N1P2) 
to  QOI.D  (NIP  I)  =  0 

1  statement  after  1001 

change  QOI.D  (Nl)  =  (ID-1)  *  QOLD  (N1M1) 
to  QOI.D  (Nl)  =  0 

1  statement  after  3000 

change  P(N1)  =  1.  -  (ID- I J  *  Q(N1MI) 
to  P(.N1)  =  1 

3  statements  after  3000 

change  Q(N1P1)  =  0  +  (2-lD)/P  (N  IP  1) 
to  Q(N1P1)  =  0 
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